KerrSchild.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 
10 #include "Options/Options.hpp"
11 #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
13 #include "Utilities/TMPL.hpp"
15 
16 /// \cond
17 namespace PUP {
18 class er; // IWYU pragma: keep
19 } // namespace PUP
20 namespace Tags {
21 template <typename Tag>
22 struct dt;
23 } // namespace Tags
24 /// \endcond
25 
26 // IWYU pragma: no_include <pup.h>
27 
28 namespace gr {
29 namespace Solutions {
30 
31 /*!
32  * \brief Kerr black hole in Kerr-Schild coordinates
33  *
34  * \details
35  * The metric is \f$g_{\mu\nu} = \eta_{\mu\nu} + 2 H l_\mu l_\nu\f$,
36  * where \f$\eta_{\mu\nu}\f$ is the Minkowski metric, \f$H\f$ is a scalar
37  * function, and \f$l_\mu\f$ is the outgoing null vector.
38  * \f$H\f$ and \f$l_\mu\f$ are known functions of the coordinates
39  * and of the mass and spin vector.
40  *
41  * The following are input file options that can be specified:
42  * - Mass (default: 1.)
43  * - Center (default: {0,0,0})
44  * - Spin (default: {0,0,0})
45  *
46  * ## Kerr-Schild Coordinates
47  *
48  *
49  * A Kerr-Schild coordinate system is defined by
50  * \f{equation}{
51  * g_{\mu\nu} \equiv \eta_{\mu\nu} + 2 H l_\mu l_\nu,
52  * \f}
53  * where \f$H\f$ is a scalar function of the coordinates, \f$\eta_{\mu\nu}\f$
54  * is the Minkowski metric, and \f$l^\mu\f$ is a null vector. Note that the
55  * form of the metric along with the nullness of \f$l^\mu\f$ allows one to
56  * raise and lower indices of \f$l^\mu\f$ using \f$\eta_{\mu\nu}\f$, and
57  * that \f$l^t l^t = l_t l_t = l^i l_i\f$.
58  * Note also that
59  * \f{equation}{
60  * g^{\mu\nu} \equiv \eta^{\mu\nu} - 2 H l^\mu l^\nu,
61  * \f}
62  * and that \f$\sqrt{-g}=1\f$.
63  * Also, \f$l_\mu\f$ is a geodesic with respect to both the physical metric
64  * and the Minkowski metric:
65  * \f{equation}{
66  * l^\mu \partial_\mu l_\nu = l^\mu\nabla_\mu l_\nu = 0.
67  * \f}
68  *
69  *
70  * The corresponding 3+1 quantities are
71  * \f{eqnarray}{
72  * g_{i j} &=& \delta_{i j} + 2 H l_i l_j,\\
73  * g^{i j} &=& \delta^{i j} - {2 H l^i l^j \over 1+2H l^t l^t},\\
74  * {\rm det} g_{i j}&=& 1+2H l^t l^t,\\
75  * \beta^i &=& - {2 H l^t l^i \over 1+2H l^t l^t},\\
76  * N &=& \left(1+2 H l^t l^t\right)^{-1/2},\quad\hbox{(lapse)}\\
77  * \alpha &=& \left(1+2 H l^t l^t\right)^{-1},
78  * \quad\hbox{(densitized lapse)}\\
79  * K_{i j} &=& - \left(1+2 H l^t l^t\right)^{1/2}
80  * \left[l_i l_j \partial_{t} H + 2 H l_{(i} \partial_{t} l_{j)}\right]
81  * \nonumber \\
82  * &&-2\left(1+2 H l^t l^t\right)^{-1/2}
83  * \left[H l^t \partial_{(i}l_{j)} + H l_{(i}\partial_{j)}l^t
84  * + l^t l_{(i}\partial_{j)} H + 2H^2 l^t l_{(i} l^k\partial_{k}l_{j)}
85  * + H l^t l_i l_j l^k \partial_{k} H\right],\\
86  * \partial_{k}g_{i j}&=& 2 l_i l_j\partial_{k} H +
87  * 4 H l_{(i} \partial_{k}l_{j)},\\
88  * \partial_{k}N &=& -\left(1+2 H l^t l^t\right)^{-3/2}
89  * \left(l^tl^t\partial_{k}H+2Hl^t\partial_{k}l^t\right),\\
90  * \partial_{k}\beta^i &=& - 2\left(1+2H l^t l^t\right)^{-1}
91  * \left(l^tl^i\partial_{k}H+Hl^t\partial_{k}l^i+Hl^i\partial_{k}l^t\right)
92  * + 4 H l^t l^i \left(1+2H l^t l^t\right)^{-2}
93  * \left(l^tl^t\partial_{k}H+2Hl^t\partial_{k}l^t\right),\\
94  * \Gamma^k{}_{i j}&=& -\delta^{k m}\left(l_i l_j \partial_{m}H
95  * + 2l_{(i} \partial_{m}l_{j)} \right)
96  * + 2 H l_{(i}\partial_{j)} l^k
97  * \nonumber \\
98  * &&+\left(1+2 H l^t l^t\right)^{-1}
99  * \left[2 l^k l_{(i}\partial_{j)} H
100  * +2 H l_i l_j l^k l^m \partial_{m}H
101  * +2 H l^k \partial_{(i}l_{j)}
102  * +4 H^2 l^k l_{(i} (l^m \partial_{m}l_{j)}
103  * -\partial_{j)} l^t)
104  * \right].
105  * \f}
106  * Note that \f$l^i\f$ is **not** equal to \f$g^{i j} l_j\f$; it is equal
107  * to \f${}^{(4)}g^{i \mu} l_\mu\f$.
108  *
109  * ## Kerr Spacetime
110  *
111  * ### Spin in the z direction
112  *
113  * Assume Cartesian coordinates \f$(t,x,y,z)\f$. Then for stationary Kerr
114  * spacetime with mass \f$M\f$ and angular momentum \f$a M\f$
115  * in the \f$z\f$ direction,
116  * \f{eqnarray}{
117  * H &=& {M r^3 \over r^4 + a^2 z^2},\\
118  * l_\mu &=&
119  * \left(1,{rx+ay\over r^2+a^2},{ry-ax\over r^2+a^2},{z\over r}\right),
120  * \f}
121  * where \f$r\f$ is defined by
122  * \f{equation}{
123  * \label{eq:rdefinition1}
124  * {x^2+y^2\over a^2+r^2} + {z^2\over r^2} = 1,
125  * \f}
126  * or equivalently,
127  * \f{equation}{
128  * r^2 = {1\over 2}(x^2 + y^2 + z^2 - a^2)
129  * + \left({1\over 4}(x^2 + y^2 + z^2 - a^2)^2 + a^2 z^2\right)^{1/2}.
130  * \f}
131  *
132  * Possibly useful formula:
133  * \f{equation}{
134  * \partial_{i} r = {x_i + z \delta_{i z} \displaystyle {a^2\over r^2} \over
135  * 2 r\left(1 - \displaystyle {x^2 + y^2 + z^2 - a^2\over 2 r^2}\right)}.
136  * \f}
137  *
138  * ### Spin in an arbitrary direction
139  *
140  * For arbitrary spin direction, let \f$\vec{x}\equiv (x,y,z)\f$ and
141  * \f$\vec{a}\f$ be a flat-space three-vector with magnitude-squared
142  * (\f$\delta_{ij}\f$ norm) equal to \f$a^2\f$.
143  * Then the Kerr-Schild quantities for Kerr spacetime are:
144  * \f{eqnarray}{
145  * H &=& {M r^3 \over r^4 + (\vec{a}\cdot\vec{x})^2},\\
146  * \vec{l} &=& {r\vec{x}-\vec{a}\times\vec{x}+(\vec{a}\cdot\vec{x})\vec{a}/r
147  * \over r^2+a^2 },\\
148  * l_t &=& 1,\\
149  * \label{eq:rdefinition2}
150  * r^2 &=& {1\over 2}(\vec{x}\cdot\vec{x}-a^2)
151  * + \left({1\over 4}(\vec{x}\cdot\vec{x}-a^2)^2
152  * + (\vec{a}\cdot\vec{x})^2\right)^{1/2},
153  * \f}
154  * where \f$\vec{l}\equiv (l_x,l_y,l_z)\f$, and
155  * all dot and cross products are evaluated as flat-space 3-vector operations.
156  *
157  * Possibly useful formulae:
158  * \f{equation}{
159  * \partial_{i} r = {x_i + (\vec{a}\cdot\vec{x})a_i/r^2 \over
160  * 2 r\left(1 - \displaystyle {\vec{x}\cdot\vec{x}-a^2\over 2 r^2}\right)},
161  * \f}
162  * \f{equation}{
163  * {\partial_{i} H \over H} =
164  * {3\partial_{i}r\over r} - {4 r^3 \partial_{i}r +
165  * 2(\vec{a}\cdot\vec{x})\vec{a}
166  * \over r^4 + (\vec{a}\cdot\vec{x})^2},
167  * \f}
168  * \f{equation}{
169  * (r^2+a^2)\partial_{j} l_i =
170  * (x_i-2 r l_i-(\vec{a}\cdot\vec{x})a_i/r^2)\partial_{j}r
171  * + r\delta_{ij} + a_i a_j/r + \epsilon^{ijk} a_k.
172  * \f}
173  *
174  * ## Cartesian and Spherical Coordinates for Kerr
175  *
176  * The Kerr-Schild coordinates are defined in terms of the Cartesian
177  * coordinates \f$(x,y,z)\f$. If one wishes to express Kerr-Schild
178  * coordinates in terms of the spherical polar coordinates
179  * \f$(\tilde{r},\theta,\phi)\f$ then one can make the obvious and
180  * usual transformation
181  * \f{equation}{
182  * \label{eq:sphertocartsimple}
183  * x=\tilde{r}\sin\theta\cos\phi,\quad
184  * y=\tilde{r}\sin\theta\sin\phi,\quad
185  * z=\tilde{r}\cos\theta.
186  * \f}
187  *
188  * This is simple, and has the advantage that in this coordinate system
189  * for \f$M\to0\f$, Kerr spacetime becomes Minkowski space in spherical
190  * coordinates \f$(\tilde{r},\theta,\phi)\f$. However, the disadvantage is
191  * that the horizon of a Kerr hole is **not** located at constant
192  * \f$\tilde{r}\f$, but is located instead at constant \f$r\f$,
193  * where \f$r\f$ is the radial
194  * Boyer-Lindquist coordinate defined in (\f$\ref{eq:rdefinition2}\f$).
195  *
196  * For spin in the \f$z\f$ direction, one could use the transformation
197  * \f{equation}{
198  * x=\sqrt{r^2+a^2}\sin\theta\cos\phi,\quad
199  * y=\sqrt{r^2+a^2}\sin\theta\sin\phi,\quad
200  * z=r\cos\theta.
201  * \f}
202  * In this case, for \f$M\to0\f$, Kerr spacetime becomes Minkowski space in
203  * spheroidal coordinates, but now the horizon is on a constant-coordinate
204  * surface.
205  *
206  * Right now we use (\f$\ref{eq:sphertocartsimple}\f$), but we may
207  * wish to use the other transformation in the future.
208  */
209 class KerrSchild {
210  public:
211  struct Mass {
212  using type = double;
213  static constexpr OptionString help = {"Mass of the black hole"};
214  static type default_value() noexcept { return 1.; }
215  static type lower_bound() noexcept { return 0.; }
216  };
217  struct Spin {
218  using type = std::array<double, 3>;
219  static constexpr OptionString help = {
220  "The [x,y,z] dimensionless spin of the black hole"};
221  static type default_value() noexcept { return {{0., 0., 0.}}; }
222  };
223  struct Center {
224  using type = std::array<double, 3>;
225  static constexpr OptionString help = {
226  "The [x,y,z] center of the black hole"};
227  static type default_value() noexcept { return {{0., 0., 0.}}; }
228  };
229  using options = tmpl::list<Mass, Spin, Center>;
230  static constexpr OptionString help{"Black hole in Kerr-Schild coordinates"};
231 
232  KerrSchild(double mass, Spin::type dimensionless_spin, Center::type center,
233  const OptionContext& context = {});
234 
235  explicit KerrSchild(CkMigrateMessage* /*unused*/) noexcept {}
236 
237  KerrSchild() = default;
238  KerrSchild(const KerrSchild& /*rhs*/) = default;
239  KerrSchild& operator=(const KerrSchild& /*rhs*/) = default;
240  KerrSchild(KerrSchild&& /*rhs*/) noexcept = default;
241  KerrSchild& operator=(KerrSchild&& /*rhs*/) noexcept = default;
242  ~KerrSchild() = default;
243 
244  template <typename DataType>
245  using DerivLapse = ::Tags::deriv<gr::Tags::Lapse<DataType>, tmpl::size_t<3>,
247  template <typename DataType>
248  using DerivShift =
249  ::Tags::deriv<gr::Tags::Shift<3, Frame::Inertial, DataType>,
250  tmpl::size_t<3>, Frame::Inertial>;
251  template <typename DataType>
252  using DerivSpatialMetric =
253  ::Tags::deriv<gr::Tags::SpatialMetric<3, Frame::Inertial, DataType>,
254  tmpl::size_t<3>, Frame::Inertial>;
255  template <typename DataType>
256  using tags = tmpl::list<
258  DerivLapse<DataType>, gr::Tags::Shift<3, Frame::Inertial, DataType>,
260  DerivShift<DataType>,
263  DerivSpatialMetric<DataType>, gr::Tags::SqrtDetSpatialMetric<DataType>,
266 
267  template <typename DataType>
268  tuples::tagged_tuple_from_typelist<tags<DataType>> variables(
269  const tnsr::I<DataType, 3>& x, double t, tags<DataType> /*meta*/) const
270  noexcept;
271 
272  // clang-tidy: no runtime references
273  void pup(PUP::er& p) noexcept; // NOLINT
274 
275  SPECTRE_ALWAYS_INLINE double mass() const noexcept { return mass_; }
276  SPECTRE_ALWAYS_INLINE const std::array<double, 3>& center() const noexcept {
277  return center_;
278  }
279  SPECTRE_ALWAYS_INLINE const std::array<double, 3>& dimensionless_spin() const
280  noexcept {
281  return dimensionless_spin_;
282  }
283 
284  private:
285  double mass_{1.0};
286  std::array<double, 3> dimensionless_spin_{{0.0, 0.0, 0.0}},
287  center_{{0.0, 0.0, 0.0}};
288 };
289 
290 SPECTRE_ALWAYS_INLINE bool operator==(const KerrSchild& lhs,
291  const KerrSchild& rhs) noexcept {
292  return lhs.mass() == rhs.mass() and
293  lhs.dimensionless_spin() == rhs.dimensionless_spin() and
294  lhs.center() == rhs.center();
295 }
296 
297 SPECTRE_ALWAYS_INLINE bool operator!=(const KerrSchild& lhs,
298  const KerrSchild& rhs) noexcept {
299  return not(lhs == rhs);
300 }
301 } // namespace Solutions
302 } // namespace gr
Definition: Strahlkorper.hpp:14
Defines class tuples::TaggedTuple.
Definition: Tags.hpp:36
Defines classes and functions for making classes creatable from input files.
Definition: Tags.hpp:46
Definition: Tags.hpp:31
Prefix indicating a time derivative.
Definition: Prefixes.hpp:28
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:26
Holds functions related to general relativity.
Definition: KerrHorizon.cpp:14
Definition: KerrSchild.hpp:217
#define SPECTRE_ALWAYS_INLINE
Always inline a function. Only use this if you benchmarked the code.
Definition: ForceInline.hpp:20
Definition: KerrSchild.hpp:223
Definition: Tags.hpp:26
Definition: KerrSchild.hpp:211
Defines functions computing partial derivatives.
Definition: Tags.hpp:41
Defines macro to always inline a function.
Definition: DataBoxTag.hpp:29
Defines a list of useful type aliases for tensors.
Information about the nested operations being performed by the parser, for use in printing errors...
Definition: Options.hpp:35
Wraps the template metaprogramming library used (brigand)
Kerr black hole in Kerr-Schild coordinates.
Definition: KerrSchild.hpp:209
Definition: IndexType.hpp:44
Definition: Tags.hpp:96