Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <limits>
8 :
9 : #include "DataStructures/TaggedTuple.hpp"
10 : #include "DataStructures/Tensor/TypeAliases.hpp"
11 : #include "Options/Context.hpp"
12 : #include "Options/String.hpp"
13 : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
14 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Solutions.hpp"
15 : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
16 : #include "Utilities/ForceInline.hpp"
17 : #include "Utilities/TMPL.hpp"
18 :
19 : /// \cond
20 : namespace PUP {
21 : class er;
22 : } // namespace PUP
23 : namespace Tags {
24 : template <typename Tag>
25 : struct dt;
26 : template <typename Tag, typename Dim, typename Frame>
27 : struct deriv;
28 : } // namespace Tags
29 : /// \endcond
30 :
31 : namespace gr::Solutions {
32 :
33 : /*!
34 : * \brief Kerr initial data in the puncture radial coordinate of Liu, Etienne,
35 : * and Shapiro \cite Liu2009.
36 : *
37 : * \details
38 : * This class represents a Kerr black hole of mass \f$M\f$ and dimensionless
39 : * spin \f$\chi\f$ (with the spin along the \f$+z\f$ axis and the puncture at
40 : * the coordinate origin) in the puncture radial coordinate \f$r\f$ introduced
41 : * in \cite Liu2009. Defining the spin \f$a = \chi M\f$ and the horizon radii
42 : * \f$r_\pm = M \pm \sqrt{M^2 - a^2}\f$ (so that \f$r_+ r_- = a^2\f$ and
43 : * \f$r_+ + r_- = 2M\f$), the Boyer-Lindquist radius \f$r_\mathrm{BL}\f$ is
44 : * related to \f$r\f$ by
45 : *
46 : * \f{align}{
47 : * r_\mathrm{BL} = r\left(1 + \frac{r_+}{4r}\right)^2
48 : * = r + \frac{r_+}{2} + \frac{r_+^2}{16 r} .
49 : * \f}
50 : *
51 : * The horizon (throat) sits at the coordinate radius \f$r = r_+/4\f$, which
52 : * stays finite (\f$\to M/4\f$) as \f$|\chi| \to 1\f$. The coordinate covers the
53 : * black hole exterior twice: the two sheets \f$r \gtrless r_+/4\f$ are joined
54 : * at the throat, and \f$r_\mathrm{BL} \ge r_+\f$ everywhere, so the interior
55 : * \f$r_\mathrm{BL} < r_+\f$ is never entered. Since the metric component
56 : * \f$\gamma_{rr}\f$ diverges at the throat like
57 : * \f$1/\sqrt{M^2 - a^2}\f$ as \f$|\chi| \to 1\f$ (the well-known infinite
58 : * proper throat of extremal Kerr), the class requires \f$|\chi| < 1\f$
59 : * strictly.
60 : *
61 : * The Boyer-Lindquist scalars are
62 : *
63 : * \f{align}{
64 : * \Sigma &= r_\mathrm{BL}^2 + \frac{a^2 z^2}{r^2}, \\
65 : * \Delta &= r_\mathrm{BL}^2 - 2 M r_\mathrm{BL} + a^2
66 : * = (r_\mathrm{BL} - r_+)(r_\mathrm{BL} - r_-), \\
67 : * A &= (r_\mathrm{BL}^2 + a^2)^2
68 : * - \Delta\, a^2 \left(1 - \frac{z^2}{r^2}\right) .
69 : * \f}
70 : *
71 : * The tensors are assembled directly in Cartesian form from the three mutually
72 : * orthogonal building blocks
73 : *
74 : * \f{align}{
75 : * n_i = \frac{x_i}{r}, \qquad
76 : * \lambda_i = (-y,\, x,\, 0), \qquad
77 : * \mu_i = (z x,\, z y,\, -(x^2 + y^2)) ,
78 : * \f}
79 : *
80 : * which are regular for \f$r > 0\f$ including on the spin axis, so no near-axis
81 : * special case is needed. The spatial metric and extrinsic curvature are
82 : *
83 : * \f{align}{
84 : * \gamma_{ij} &= c_\delta\, \delta_{ij} + c_n\, n_i n_j
85 : * + c_\lambda\, \lambda_i \lambda_j, \\
86 : * K_{ij} &= c_{n\lambda}\, (n_i \lambda_j + n_j \lambda_i)
87 : * + c_{\mu\lambda}\, (\mu_i \lambda_j + \mu_j \lambda_i),
88 : * \f}
89 : *
90 : * with the coefficient functions
91 : *
92 : * \f{align}{
93 : * c_\delta &= \frac{\Sigma}{r^2}, \qquad
94 : * c_n = \frac{\Sigma\, r_-}{r^2 (r_\mathrm{BL} - r_-)}, \qquad
95 : * c_\lambda = \frac{a^2 (\Sigma + 2 M r_\mathrm{BL})}{\Sigma\, r^4}, \\
96 : * c_{n\lambda} &= \frac{M a\, G\, \sqrt{r_\mathrm{BL}}}
97 : * {\Sigma \sqrt{A \Sigma}\; r^3 \sqrt{r_\mathrm{BL} - r_-}}, \qquad
98 : * c_{\mu\lambda} = -\frac{2 a^3 M r_\mathrm{BL}\, z\, (r - r_+/4)}
99 : * {\Sigma \sqrt{A \Sigma}\; r^6}\sqrt{\frac{r_\mathrm{BL} - r_-}{r}},
100 : * \f}
101 : *
102 : * where
103 : *
104 : * \f{align}{
105 : * G = 3 r_\mathrm{BL}^4 + 2 a^2 r_\mathrm{BL}^2 - a^4
106 : * - a^2 (r_\mathrm{BL}^2 - a^2)\left(1 - \frac{z^2}{r^2}\right) .
107 : * \f}
108 : *
109 : * The stationary shift of \cite Liu2009 Eq. (7) is purely azimuthal,
110 : * \f$\beta^\phi = -2 M a r_\mathrm{BL}/A\f$, which in Cartesian form is
111 : *
112 : * \f{align}{
113 : * \beta^i = c_\beta\, \lambda^i, \qquad
114 : * c_\beta = -\frac{2 M a\, r_\mathrm{BL}}{A} .
115 : * \f}
116 : *
117 : * The analytic lapse
118 : *
119 : * \f{align}{
120 : * \alpha = \left(r - \frac{r_+}{4}\right) g(r, z) ,
121 : * \f}
122 : *
123 : * is negative on the inner sheet \f$r < r_+/4\f$. With this lapse the
124 : * returned representation is exactly stationary on both sheets, and all time
125 : * derivative tags vanish identically. Consumers that require an
126 : * everywhere-nonnegative initial lapse (e.g. moving-puncture evolutions)
127 : * must take the absolute value, recovering \cite Liu2009 Eq. (6).
128 : *
129 : * The inverse spatial metric and \f$\sqrt{\det\gamma}\f$ are returned from
130 : * their closed forms
131 : *
132 : * \f{align}{
133 : * \gamma^{ij} &= \frac{1}{c_\delta}\, \delta^{ij}
134 : * - \frac{c_n}{c_\delta (c_\delta + c_n)}\, n^i n^j
135 : * - \frac{c_\lambda}{c_\delta (c_\delta + \varpi^2 c_\lambda)}\,
136 : * \lambda^i \lambda^j, \\
137 : * \sqrt{\det\gamma} &= \frac{1}{r^3}
138 : * \sqrt{\frac{\Sigma\, r_\mathrm{BL}\, A}{r_\mathrm{BL} - r_-}},
139 : * \f}
140 : *
141 : * with \f$\varpi^2 = x^2 + y^2\f$, rather than by numerical inversion.
142 : *
143 : * Like gr::Solutions::TrumpetSchwarzschild, quantities diverge at the
144 : * puncture point \f$r = 0\f$ itself (the second asymptotically flat end); no
145 : * clamping is applied, so the origin must not coincide with a grid point.
146 : *
147 : * This solution reduces to the standard Schwarzschild solution in
148 : * isotropic coordinates when \f$\chi = 0\f$.
149 : *
150 : * The following input file options can be specified:
151 : * - Mass (\f$M > 0\f$)
152 : * - DimensionlessSpin (\f$\chi\f$, with \f$|\chi| < 1\f$ strictly)
153 : */
154 1 : class HighSpinKerrPuncture : public MarkAsAnalyticSolution,
155 : public AnalyticSolution<3_st> {
156 : private:
157 : template <typename DataType>
158 : struct IntermediateVars;
159 :
160 : public:
161 0 : static constexpr size_t volume_dim = 3;
162 :
163 0 : struct Mass {
164 0 : using type = double;
165 0 : static constexpr Options::String help = {"Mass of the Kerr black hole"};
166 0 : static type lower_bound() { return 0.; };
167 : };
168 :
169 0 : struct DimensionlessSpin {
170 0 : using type = double;
171 0 : static constexpr Options::String help = {
172 : "Dimensionless spin chi of the Kerr black hole, along the +z axis. "
173 : "Must satisfy |chi| < 1 strictly, since the metric degenerates at the "
174 : "throat in the extremal limit."};
175 : };
176 :
177 0 : using options = tmpl::list<Mass, DimensionlessSpin>;
178 0 : static constexpr Options::String help{
179 : "Kerr solution in the puncture radial coordinate of Liu, Etienne, and "
180 : "Shapiro (2009)."};
181 :
182 0 : HighSpinKerrPuncture(double mass, double dimensionless_spin,
183 : const Options::Context& context = {});
184 :
185 0 : HighSpinKerrPuncture() = default;
186 0 : HighSpinKerrPuncture(const HighSpinKerrPuncture& /*rhs*/) = default;
187 0 : HighSpinKerrPuncture& operator=(const HighSpinKerrPuncture& /*rhs*/) =
188 : default;
189 0 : HighSpinKerrPuncture(HighSpinKerrPuncture&& /*rhs*/) = default;
190 0 : HighSpinKerrPuncture& operator=(HighSpinKerrPuncture&& /*rhs*/) = default;
191 0 : ~HighSpinKerrPuncture() = default;
192 :
193 0 : explicit HighSpinKerrPuncture(CkMigrateMessage* /*msg*/);
194 :
195 : template <typename DataType>
196 0 : using DerivLapse = ::Tags::deriv<gr::Tags::Lapse<DataType>,
197 : tmpl::size_t<volume_dim>, Frame::Inertial>;
198 : template <typename DataType>
199 0 : using DerivShift = ::Tags::deriv<gr::Tags::Shift<DataType, volume_dim>,
200 : tmpl::size_t<volume_dim>, Frame::Inertial>;
201 : template <typename DataType>
202 0 : using DerivSpatialMetric =
203 : ::Tags::deriv<gr::Tags::SpatialMetric<DataType, volume_dim>,
204 : tmpl::size_t<volume_dim>, Frame::Inertial>;
205 :
206 : template <typename DataType>
207 0 : using tags = tmpl::list<
208 : gr::Tags::Lapse<DataType>, ::Tags::dt<gr::Tags::Lapse<DataType>>,
209 : DerivLapse<DataType>, gr::Tags::Shift<DataType, volume_dim>,
210 : ::Tags::dt<gr::Tags::Shift<DataType, volume_dim>>, DerivShift<DataType>,
211 : gr::Tags::SpatialMetric<DataType, volume_dim>,
212 : ::Tags::dt<gr::Tags::SpatialMetric<DataType, volume_dim>>,
213 : DerivSpatialMetric<DataType>, gr::Tags::SqrtDetSpatialMetric<DataType>,
214 : gr::Tags::ExtrinsicCurvature<DataType, volume_dim>,
215 : gr::Tags::InverseSpatialMetric<DataType, volume_dim>>;
216 :
217 : template <typename DataType, typename... Tags>
218 0 : tuples::TaggedTuple<Tags...> variables(
219 : const tnsr::I<DataType, volume_dim, Frame::Inertial>& x, double t,
220 : tmpl::list<Tags...> /*meta*/) const {
221 : const auto& vars =
222 : IntermediateVars<DataType>{mass_, dimensionless_spin_, x};
223 : return {get<Tags>(variables(x, t, vars, tmpl::list<Tags>{}))...};
224 : }
225 :
226 : // NOLINTNEXTLINE(google-runtime-references)
227 0 : void pup(PUP::er& p);
228 :
229 0 : SPECTRE_ALWAYS_INLINE double mass() const { return mass_; }
230 0 : SPECTRE_ALWAYS_INLINE double dimensionless_spin() const {
231 : return dimensionless_spin_;
232 : }
233 :
234 : private:
235 : template <typename DataType>
236 0 : auto variables(const tnsr::I<DataType, volume_dim, Frame::Inertial>& x,
237 : double t, const IntermediateVars<DataType>& vars,
238 : tmpl::list<gr::Tags::Lapse<DataType>> /*meta*/) const
239 : -> tuples::TaggedTuple<gr::Tags::Lapse<DataType>>;
240 :
241 : template <typename DataType>
242 0 : auto variables(const tnsr::I<DataType, volume_dim, Frame::Inertial>& x,
243 : double t, const IntermediateVars<DataType>& vars,
244 : tmpl::list<::Tags::dt<gr::Tags::Lapse<DataType>>> /*meta*/)
245 : const -> tuples::TaggedTuple<::Tags::dt<gr::Tags::Lapse<DataType>>>;
246 :
247 : template <typename DataType>
248 0 : auto variables(const tnsr::I<DataType, volume_dim, Frame::Inertial>& x,
249 : double t, const IntermediateVars<DataType>& vars,
250 : tmpl::list<DerivLapse<DataType>> /*meta*/) const
251 : -> tuples::TaggedTuple<DerivLapse<DataType>>;
252 :
253 : template <typename DataType>
254 0 : auto variables(const tnsr::I<DataType, volume_dim, Frame::Inertial>& x,
255 : double t, const IntermediateVars<DataType>& vars,
256 : tmpl::list<gr::Tags::Shift<DataType, volume_dim>> /*meta*/)
257 : const -> tuples::TaggedTuple<gr::Tags::Shift<DataType, volume_dim>>;
258 :
259 : template <typename DataType>
260 0 : auto variables(
261 : const tnsr::I<DataType, volume_dim, Frame::Inertial>& x, double t,
262 : const IntermediateVars<DataType>& vars,
263 : tmpl::list<::Tags::dt<gr::Tags::Shift<DataType, volume_dim>>> /*meta*/)
264 : const
265 : -> tuples::TaggedTuple<::Tags::dt<gr::Tags::Shift<DataType, volume_dim>>>;
266 :
267 : template <typename DataType>
268 0 : auto variables(const tnsr::I<DataType, volume_dim, Frame::Inertial>& x,
269 : double t, const IntermediateVars<DataType>& vars,
270 : tmpl::list<DerivShift<DataType>> /*meta*/) const
271 : -> tuples::TaggedTuple<DerivShift<DataType>>;
272 :
273 : template <typename DataType>
274 0 : auto variables(
275 : const tnsr::I<DataType, volume_dim, Frame::Inertial>& x, double t,
276 : const IntermediateVars<DataType>& vars,
277 : tmpl::list<gr::Tags::SpatialMetric<DataType, volume_dim>> /*meta*/) const
278 : -> tuples::TaggedTuple<gr::Tags::SpatialMetric<DataType, volume_dim>>;
279 :
280 : template <typename DataType>
281 0 : auto variables(
282 : const tnsr::I<DataType, volume_dim, Frame::Inertial>& x, double t,
283 : const IntermediateVars<DataType>& vars,
284 : tmpl::list<
285 : ::Tags::dt<gr::Tags::SpatialMetric<DataType, volume_dim>>> /*meta*/)
286 : const -> tuples::TaggedTuple<
287 : ::Tags::dt<gr::Tags::SpatialMetric<DataType, volume_dim>>>;
288 :
289 : template <typename DataType>
290 0 : auto variables(const tnsr::I<DataType, volume_dim, Frame::Inertial>& x,
291 : double t, const IntermediateVars<DataType>& vars,
292 : tmpl::list<DerivSpatialMetric<DataType>> /*meta*/) const
293 : -> tuples::TaggedTuple<DerivSpatialMetric<DataType>>;
294 :
295 : template <typename DataType>
296 0 : auto variables(const tnsr::I<DataType, volume_dim, Frame::Inertial>& x,
297 : double t, const IntermediateVars<DataType>& vars,
298 : tmpl::list<gr::Tags::SqrtDetSpatialMetric<DataType>> /*meta*/)
299 : const -> tuples::TaggedTuple<gr::Tags::SqrtDetSpatialMetric<DataType>>;
300 :
301 : template <typename DataType>
302 0 : auto variables(
303 : const tnsr::I<DataType, volume_dim, Frame::Inertial>& x, double t,
304 : const IntermediateVars<DataType>& vars,
305 : tmpl::list<gr::Tags::ExtrinsicCurvature<DataType, volume_dim>> /*meta*/)
306 : const -> tuples::TaggedTuple<
307 : gr::Tags::ExtrinsicCurvature<DataType, volume_dim>>;
308 :
309 : template <typename DataType>
310 0 : auto variables(
311 : const tnsr::I<DataType, volume_dim, Frame::Inertial>& x, double t,
312 : const IntermediateVars<DataType>& vars,
313 : tmpl::list<gr::Tags::InverseSpatialMetric<DataType, volume_dim>> /*meta*/)
314 : const -> tuples::TaggedTuple<
315 : gr::Tags::InverseSpatialMetric<DataType, volume_dim>>;
316 :
317 : // Intermediate quantities, computed once per call to variables().
318 : // Construct the radial building blocks r and z, the derived
319 : // constants r_+, r_-, the Boyer-Lindquist radius r_BL and its r-derivative,
320 : // the Boyer-Lindquist scalars Sigma, Delta, A, G and their (r, z) partials,
321 : // then the coefficient functions of the Cartesian assembly and their (r, z)
322 : // partials.
323 : template <typename DataType>
324 0 : struct IntermediateVars {
325 0 : IntermediateVars(double mass, double dimensionless_spin,
326 : const tnsr::I<DataType, volume_dim, Frame::Inertial>& x);
327 :
328 : // Coordinate building blocks
329 0 : DataType r{};
330 0 : DataType z{};
331 0 : DataType one_over_r{};
332 :
333 : // Boyer-Lindquist radius and its r-derivative
334 0 : DataType r_bl{};
335 0 : DataType d_r_bl_d_r{};
336 :
337 : // Boyer-Lindquist scalars
338 0 : DataType sigma{};
339 0 : DataType delta{};
340 0 : DataType a_capital{};
341 0 : DataType g_capital{};
342 0 : DataType r_bl_minus_r_minus{};
343 0 : DataType r_minus_r_plus_over_four{};
344 :
345 : // (r, z) partials of the Boyer-Lindquist scalars
346 0 : DataType d_sigma_d_r{};
347 0 : DataType d_sigma_d_z{};
348 0 : DataType d_a_capital_d_r{};
349 0 : DataType d_a_capital_d_z{};
350 0 : DataType d_g_capital_d_r{};
351 0 : DataType d_g_capital_d_z{};
352 :
353 : // Cartesian coefficient functions and their (r, z) partials
354 0 : DataType c_delta{};
355 0 : DataType c_n{};
356 0 : DataType c_lambda{};
357 0 : DataType c_n_lambda{};
358 : // c_mu_lambda = (r - r_+/4) * z * c_mu_lambda_base. The two vanishing
359 : // factors (r - r_+/4) at the throat and z on the equator are kept
360 : // explicit so the coefficient and its derivatives avoid 0/0.
361 0 : DataType c_mu_lambda_base{};
362 0 : DataType c_mu_lambda{};
363 0 : DataType g_lapse{};
364 0 : DataType c_beta{};
365 :
366 0 : DataType d_c_delta_d_r{};
367 0 : DataType d_c_delta_d_z{};
368 0 : DataType d_c_n_d_r{};
369 0 : DataType d_c_n_d_z{};
370 0 : DataType d_c_lambda_d_r{};
371 0 : DataType d_c_lambda_d_z{};
372 0 : DataType d_c_n_lambda_d_r{};
373 0 : DataType d_c_n_lambda_d_z{};
374 0 : DataType d_c_mu_lambda_d_r{};
375 0 : DataType d_c_mu_lambda_d_z{};
376 0 : DataType d_g_lapse_d_r{};
377 0 : DataType d_g_lapse_d_z{};
378 0 : DataType d_c_beta_d_r{};
379 0 : DataType d_c_beta_d_z{};
380 :
381 : // Derived scalar constants of the black hole
382 0 : double mass;
383 0 : double spin_a;
384 0 : double r_plus;
385 0 : double r_minus;
386 : };
387 :
388 0 : double mass_{std::numeric_limits<double>::signaling_NaN()};
389 0 : double dimensionless_spin_{std::numeric_limits<double>::signaling_NaN()};
390 : };
391 :
392 0 : bool operator==(const HighSpinKerrPuncture& lhs,
393 : const HighSpinKerrPuncture& rhs);
394 0 : bool operator!=(const HighSpinKerrPuncture& lhs,
395 : const HighSpinKerrPuncture& rhs);
396 : } // namespace gr::Solutions
|
