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 <pup.h>
9 :
10 : #include "DataStructures/CachedTempBuffer.hpp"
11 : #include "DataStructures/Tags/TempTensor.hpp"
12 : #include "DataStructures/Tensor/TypeAliases.hpp"
13 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
14 : #include "Options/Context.hpp"
15 : #include "Options/String.hpp"
16 : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
17 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Solutions.hpp"
18 : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
19 : #include "Utilities/ContainerHelpers.hpp"
20 : #include "Utilities/ForceInline.hpp"
21 : #include "Utilities/MakeArray.hpp"
22 : #include "Utilities/TMPL.hpp"
23 : #include "Utilities/TaggedTuple.hpp"
24 :
25 : /// \cond
26 : namespace Tags {
27 : template <typename Tag>
28 : struct dt;
29 : } // namespace Tags
30 : namespace gsl {
31 : template <class T>
32 : class not_null;
33 : } // namespace gsl
34 : /// \endcond
35 :
36 : namespace gr {
37 : namespace Solutions {
38 :
39 : /*!
40 : * \brief Kerr black hole in Kerr-Schild coordinates
41 : *
42 : * \details
43 : * The metric is \f$g_{\mu\nu} = \eta_{\mu\nu} + 2 H l_\mu l_\nu\f$,
44 : * where \f$\eta_{\mu\nu}\f$ is the Minkowski metric, \f$H\f$ is a scalar
45 : * function, and \f$l_\mu\f$ is the outgoing null vector.
46 : * \f$H\f$ and \f$l_\mu\f$ are known functions of the coordinates
47 : * and of the mass and spin vector.
48 : *
49 : * The following are input file options that can be specified:
50 : * - Mass (default: 1.)
51 : * - Center (default: {0,0,0})
52 : * - Spin (default: {0,0,0})
53 : *
54 : * ## Kerr-Schild Coordinates
55 : *
56 : *
57 : * A Kerr-Schild coordinate system is defined by
58 : * \f{equation}{
59 : * g_{\mu\nu} \equiv \eta_{\mu\nu} + 2 H l_\mu l_\nu,
60 : * \f}
61 : * where \f$H\f$ is a scalar function of the coordinates, \f$\eta_{\mu\nu}\f$
62 : * is the Minkowski metric, and \f$l^\mu\f$ is a null vector. Note that the
63 : * form of the metric along with the nullness of \f$l^\mu\f$ allows one to
64 : * raise and lower indices of \f$l^\mu\f$ using \f$\eta_{\mu\nu}\f$, and
65 : * that \f$l^t l^t = l_t l_t = l^i l_i\f$.
66 : * Note also that
67 : * \f{equation}{
68 : * g^{\mu\nu} \equiv \eta^{\mu\nu} - 2 H l^\mu l^\nu,
69 : * \f}
70 : * and that \f$\sqrt{-g}=1\f$.
71 : * Also, \f$l_\mu\f$ is a geodesic with respect to both the physical metric
72 : * and the Minkowski metric:
73 : * \f{equation}{
74 : * l^\mu \partial_\mu l_\nu = l^\mu\nabla_\mu l_\nu = 0.
75 : * \f}
76 : *
77 : *
78 : * The corresponding 3+1 quantities are
79 : * \f{eqnarray}{
80 : * g_{i j} &=& \delta_{i j} + 2 H l_i l_j,\\
81 : * g^{i j} &=& \delta^{i j} - {2 H l^i l^j \over 1+2H l^t l^t},\\
82 : * {\rm det} g_{i j}&=& 1+2H l^t l^t,\\
83 : * \partial_k ({\rm det} g_{i j})&=& 2 l^t l^t \partial_k H,\\
84 : * \beta^i &=& - {2 H l^t l^i \over 1+2H l^t l^t},\\
85 : * N &=& \left(1+2 H l^t l^t\right)^{-1/2},\quad\hbox{(lapse)}\\
86 : * \alpha &=& \left(1+2 H l^t l^t\right)^{-1},
87 : * \quad\hbox{(densitized lapse)}\\
88 : * K_{i j} &=& - \left(1+2 H l^t l^t\right)^{1/2}
89 : * \left[l_i l_j \partial_{t} H + 2 H l_{(i} \partial_{t} l_{j)}\right]
90 : * \nonumber \\
91 : * &&-2\left(1+2 H l^t l^t\right)^{-1/2}
92 : * \left[H l^t \partial_{(i}l_{j)} + H l_{(i}\partial_{j)}l^t
93 : * + l^t l_{(i}\partial_{j)} H + 2H^2 l^t l_{(i} l^k\partial_{k}l_{j)}
94 : * + H l^t l_i l_j l^k \partial_{k} H\right],\\
95 : * \partial_{k}g_{i j}&=& 2 l_i l_j\partial_{k} H +
96 : * 4 H l_{(i} \partial_{k}l_{j)},\\
97 : * \partial_{k}N &=& -\left(1+2 H l^t l^t\right)^{-3/2}
98 : * \left(l^tl^t\partial_{k}H+2Hl^t\partial_{k}l^t\right),\\
99 : * \partial_{k}\beta^i &=& - 2\left(1+2H l^t l^t\right)^{-1}
100 : * \left(l^tl^i\partial_{k}H+Hl^t\partial_{k}l^i+Hl^i\partial_{k}l^t\right)
101 : * + 4 H l^t l^i \left(1+2H l^t l^t\right)^{-2}
102 : * \left(l^tl^t\partial_{k}H+2Hl^t\partial_{k}l^t\right),\\
103 : * \Gamma^k{}_{i j}&=& -\delta^{k m}\left(l_i l_j \partial_{m}H
104 : * + 2l_{(i} \partial_{m}l_{j)} \right)
105 : * + 2 H l_{(i}\partial_{j)} l^k
106 : * \nonumber \\
107 : * &&+\left(1+2 H l^t l^t\right)^{-1}
108 : * \left[2 l^k l_{(i}\partial_{j)} H
109 : * +2 H l_i l_j l^k l^m \partial_{m}H
110 : * +2 H l^k \partial_{(i}l_{j)}
111 : * +4 H^2 l^k l_{(i} (l^m \partial_{m}l_{j)}
112 : * -\partial_{j)} l^t)
113 : * \right].
114 : * \f}
115 : * Note that \f$l^i\f$ is **not** equal to \f$g^{i j} l_j\f$; it is equal
116 : * to \f${}^{(4)}g^{i \mu} l_\mu\f$.
117 : *
118 : * ## Kerr Spacetime
119 : *
120 : * ### Spin in the z direction
121 : *
122 : * Assume Cartesian coordinates \f$(t,x,y,z)\f$. Then for stationary Kerr
123 : * spacetime with mass \f$M\f$ and angular momentum \f$a M\f$
124 : * in the \f$z\f$ direction,
125 : * \f{eqnarray}{
126 : * H &=& {M r^3 \over r^4 + a^2 z^2},\\
127 : * l_\mu &=&
128 : * \left(1,{rx+ay\over r^2+a^2},{ry-ax\over r^2+a^2},{z\over r}\right),
129 : * \f}
130 : * where \f$r\f$ is defined by
131 : * \f{equation}{
132 : * \label{eq:rdefinition1}
133 : * {x^2+y^2\over a^2+r^2} + {z^2\over r^2} = 1,
134 : * \f}
135 : * or equivalently,
136 : * \f{equation}{
137 : * r^2 = {1\over 2}(x^2 + y^2 + z^2 - a^2)
138 : * + \left({1\over 4}(x^2 + y^2 + z^2 - a^2)^2 + a^2 z^2\right)^{1/2}.
139 : * \f}
140 : *
141 : * Possibly useful formula:
142 : * \f{equation}{
143 : * \partial_{i} r = {x_i + z \delta_{i z} \displaystyle {a^2\over r^2} \over
144 : * 2 r\left(1 - \displaystyle {x^2 + y^2 + z^2 - a^2\over 2 r^2}\right)}.
145 : * \f}
146 : *
147 : * ### Spin in an arbitrary direction
148 : *
149 : * For arbitrary spin direction, let \f$\vec{x}\equiv (x,y,z)\f$ and
150 : * \f$\vec{a}\f$ be a flat-space three-vector with magnitude-squared
151 : * (\f$\delta_{ij}\f$ norm) equal to \f$a^2\f$.
152 : * Then the Kerr-Schild quantities for Kerr spacetime are:
153 : * \f{eqnarray}{
154 : * H &=& {M r^3 \over r^4 + (\vec{a}\cdot\vec{x})^2},\\
155 : * \vec{l} &=& {r\vec{x}-\vec{a}\times\vec{x}+(\vec{a}\cdot\vec{x})\vec{a}/r
156 : * \over r^2+a^2 },\\
157 : * l_t &=& 1,\\
158 : * \label{eq:rdefinition2}
159 : * r^2 &=& {1\over 2}(\vec{x}\cdot\vec{x}-a^2)
160 : * + \left({1\over 4}(\vec{x}\cdot\vec{x}-a^2)^2
161 : * + (\vec{a}\cdot\vec{x})^2\right)^{1/2},
162 : * \f}
163 : * where \f$\vec{l}\equiv (l_x,l_y,l_z)\f$, and
164 : * all dot and cross products are evaluated as flat-space 3-vector operations.
165 : *
166 : * Possibly useful formulae:
167 : * \f{equation}{
168 : * \partial_{i} r = {x_i + (\vec{a}\cdot\vec{x})a_i/r^2 \over
169 : * 2 r\left(1 - \displaystyle {\vec{x}\cdot\vec{x}-a^2\over 2 r^2}\right)},
170 : * \f}
171 : * \f{equation}{
172 : * {\partial_{i} H \over H} =
173 : * {3\partial_{i}r\over r} - {4 r^3 \partial_{i}r +
174 : * 2(\vec{a}\cdot\vec{x})\vec{a}
175 : * \over r^4 + (\vec{a}\cdot\vec{x})^2},
176 : * \f}
177 : * \f{equation}{
178 : * (r^2+a^2)\partial_{j} l_i =
179 : * (x_i-2 r l_i-(\vec{a}\cdot\vec{x})a_i/r^2)\partial_{j}r
180 : * + r\delta_{ij} + a_i a_j/r + \epsilon^{ijk} a_k.
181 : * \f}
182 : *
183 : * ## Cartesian and Spherical Coordinates for Kerr
184 : *
185 : * The Kerr-Schild coordinates are defined in terms of the Cartesian
186 : * coordinates \f$(x,y,z)\f$. If one wishes to express Kerr-Schild
187 : * coordinates in terms of the spherical polar coordinates
188 : * \f$(\tilde{r},\theta,\phi)\f$ then one can make the obvious and
189 : * usual transformation
190 : * \f{equation}{
191 : * \label{eq:sphertocartsimple}
192 : * x=\tilde{r}\sin\theta\cos\phi,\quad
193 : * y=\tilde{r}\sin\theta\sin\phi,\quad
194 : * z=\tilde{r}\cos\theta.
195 : * \f}
196 : *
197 : * This is simple, and has the advantage that in this coordinate system
198 : * for \f$M\to0\f$, Kerr spacetime becomes Minkowski space in spherical
199 : * coordinates \f$(\tilde{r},\theta,\phi)\f$. However, the disadvantage is
200 : * that the horizon of a Kerr hole is **not** located at constant
201 : * \f$\tilde{r}\f$, but is located instead at constant \f$r\f$,
202 : * where \f$r\f$ is the radial
203 : * Boyer-Lindquist coordinate defined in (\f$\ref{eq:rdefinition2}\f$).
204 : *
205 : * For spin in the \f$z\f$ direction, one could use the transformation
206 : * \f{equation}{
207 : * x=\sqrt{r^2+a^2}\sin\theta\cos\phi,\quad
208 : * y=\sqrt{r^2+a^2}\sin\theta\sin\phi,\quad
209 : * z=r\cos\theta.
210 : * \f}
211 : * In this case, for \f$M\to0\f$, Kerr spacetime becomes Minkowski space in
212 : * spheroidal coordinates, but now the horizon is on a constant-coordinate
213 : * surface.
214 : *
215 : * Right now we use (\f$\ref{eq:sphertocartsimple}\f$), but we may
216 : * wish to use the other transformation in the future.
217 : */
218 1 : class KerrSchild : public AnalyticSolution<3_st>,
219 : public MarkAsAnalyticSolution {
220 : public:
221 0 : struct Mass {
222 0 : using type = double;
223 0 : static constexpr Options::String help = {"Mass of the black hole"};
224 0 : static type lower_bound() { return 0.; }
225 : };
226 0 : struct Spin {
227 0 : using type = std::array<double, volume_dim>;
228 0 : static constexpr Options::String help = {
229 : "The [x,y,z] dimensionless spin of the black hole"};
230 : };
231 0 : struct Center {
232 0 : using type = std::array<double, volume_dim>;
233 0 : static constexpr Options::String help = {
234 : "The [x,y,z] center of the black hole"};
235 : };
236 0 : using options = tmpl::list<Mass, Spin, Center>;
237 0 : static constexpr Options::String help{
238 : "Black hole in Kerr-Schild coordinates"};
239 :
240 0 : KerrSchild(double mass, const std::array<double, 3>& dimensionless_spin,
241 : const std::array<double, 3>& center,
242 : const Options::Context& context = {});
243 :
244 0 : explicit KerrSchild(CkMigrateMessage* /*msg*/);
245 :
246 : template <typename DataType, typename Frame = Frame::Inertial>
247 0 : using tags = tmpl::flatten<tmpl::list<
248 : AnalyticSolution<3_st>::tags<DataType, Frame>,
249 : gr::Tags::DerivDetSpatialMetric<DataType, 3, Frame>,
250 : gr::Tags::TraceExtrinsicCurvature<DataType>,
251 : gr::Tags::SpatialChristoffelFirstKind<DataType, 3, Frame>,
252 : gr::Tags::SpatialChristoffelSecondKind<DataType, 3, Frame>,
253 : gr::Tags::TraceSpatialChristoffelSecondKind<DataType, 3, Frame>,
254 : gr::Tags::SpacetimeChristoffelSecondKind<DataType, 3, Frame>>>;
255 :
256 0 : KerrSchild() = default;
257 0 : KerrSchild(const KerrSchild& /*rhs*/) = default;
258 0 : KerrSchild& operator=(const KerrSchild& /*rhs*/) = default;
259 0 : KerrSchild(KerrSchild&& /*rhs*/) = default;
260 0 : KerrSchild& operator=(KerrSchild&& /*rhs*/) = default;
261 0 : ~KerrSchild() = default;
262 :
263 : template <typename DataType, typename Frame, typename... Tags>
264 0 : tuples::TaggedTuple<Tags...> variables(
265 : const tnsr::I<DataType, volume_dim, Frame>& x, double /*t*/,
266 : tmpl::list<Tags...> /*meta*/) const {
267 : static_assert(
268 : tmpl2::flat_all_v<
269 : tmpl::list_contains_v<tags<DataType, Frame>, Tags>...>,
270 : "At least one of the requested tags is not supported. The requested "
271 : "tags are listed as template parameters of the `variables` function.");
272 : IntermediateVars<DataType, Frame> cache(get_size(*x.begin()));
273 : IntermediateComputer<DataType, Frame> computer(*this, x);
274 : return {cache.get_var(computer, Tags{})...};
275 : }
276 :
277 : // NOLINTNEXTLINE(google-runtime-references)
278 0 : void pup(PUP::er& p);
279 :
280 0 : double mass() const { return mass_; }
281 0 : const std::array<double, volume_dim>& center() const { return center_; }
282 0 : const std::array<double, volume_dim>& dimensionless_spin() const {
283 : return dimensionless_spin_;
284 : }
285 0 : bool zero_spin() const { return zero_spin_; }
286 :
287 0 : struct internal_tags {
288 : template <typename DataType, typename Frame = ::Frame::Inertial>
289 0 : using x_minus_center = ::Tags::TempI<0, 3, Frame, DataType>;
290 : template <typename DataType>
291 0 : using a_dot_x = ::Tags::TempScalar<1, DataType>;
292 : template <typename DataType>
293 0 : using a_dot_x_squared = ::Tags::TempScalar<2, DataType>;
294 : template <typename DataType>
295 0 : using half_xsq_minus_asq = ::Tags::TempScalar<3, DataType>;
296 : template <typename DataType>
297 0 : using r_squared = ::Tags::TempScalar<4, DataType>;
298 : template <typename DataType>
299 0 : using r = ::Tags::TempScalar<5, DataType>;
300 : template <typename DataType>
301 0 : using a_dot_x_over_rsquared = ::Tags::TempScalar<6, DataType>;
302 : template <typename DataType>
303 0 : using deriv_log_r_denom = ::Tags::TempScalar<7, DataType>;
304 : template <typename DataType, typename Frame = ::Frame::Inertial>
305 0 : using deriv_log_r = ::Tags::Tempi<8, 3, Frame, DataType>;
306 : template <typename DataType>
307 0 : using H_denom = ::Tags::TempScalar<9, DataType>;
308 : template <typename DataType>
309 0 : using H = ::Tags::TempScalar<10, DataType>;
310 : template <typename DataType>
311 0 : using deriv_H_temp1 = ::Tags::TempScalar<11, DataType>;
312 : template <typename DataType>
313 0 : using deriv_H_temp2 = ::Tags::TempScalar<12, DataType>;
314 : template <typename DataType, typename Frame = ::Frame::Inertial>
315 0 : using deriv_H = ::Tags::Tempi<13, 3, Frame, DataType>;
316 : template <typename DataType>
317 0 : using denom = ::Tags::TempScalar<14, DataType>;
318 : template <typename DataType>
319 0 : using a_dot_x_over_r = ::Tags::TempScalar<15, DataType>;
320 : template <typename DataType, typename Frame = ::Frame::Inertial>
321 0 : using null_form = ::Tags::Tempi<16, 3, Frame, DataType>;
322 : template <typename DataType, typename Frame = ::Frame::Inertial>
323 0 : using deriv_null_form = ::Tags::Tempij<17, 3, Frame, DataType>;
324 : template <typename DataType>
325 0 : using lapse_squared = ::Tags::TempScalar<18, DataType>;
326 : template <typename DataType>
327 0 : using deriv_lapse_multiplier = ::Tags::TempScalar<19, DataType>;
328 : template <typename DataType>
329 0 : using shift_multiplier = ::Tags::TempScalar<20, DataType>;
330 : };
331 :
332 : template <typename DataType, typename Frame = ::Frame::Inertial>
333 0 : using CachedBuffer = CachedTempBuffer<
334 : internal_tags::x_minus_center<DataType, Frame>,
335 : internal_tags::a_dot_x<DataType>,
336 : internal_tags::a_dot_x_squared<DataType>,
337 : internal_tags::half_xsq_minus_asq<DataType>,
338 : internal_tags::r_squared<DataType>, internal_tags::r<DataType>,
339 : internal_tags::a_dot_x_over_rsquared<DataType>,
340 : internal_tags::deriv_log_r_denom<DataType>,
341 : internal_tags::deriv_log_r<DataType, Frame>,
342 : internal_tags::H_denom<DataType>, internal_tags::H<DataType>,
343 : internal_tags::deriv_H_temp1<DataType>,
344 : internal_tags::deriv_H_temp2<DataType>,
345 : internal_tags::deriv_H<DataType, Frame>, internal_tags::denom<DataType>,
346 : internal_tags::a_dot_x_over_r<DataType>,
347 : internal_tags::null_form<DataType, Frame>,
348 : internal_tags::deriv_null_form<DataType, Frame>,
349 : internal_tags::lapse_squared<DataType>, gr::Tags::Lapse<DataType>,
350 : internal_tags::deriv_lapse_multiplier<DataType>,
351 : internal_tags::shift_multiplier<DataType>,
352 : gr::Tags::Shift<DataType, 3, Frame>, DerivShift<DataType, Frame>,
353 : gr::Tags::SpatialMetric<DataType, 3, Frame>,
354 : gr::Tags::InverseSpatialMetric<DataType, 3, Frame>,
355 : DerivSpatialMetric<DataType, Frame>,
356 : ::Tags::dt<gr::Tags::SpatialMetric<DataType, 3, Frame>>,
357 : gr::Tags::ExtrinsicCurvature<DataType, 3, Frame>,
358 : gr::Tags::SpatialChristoffelFirstKind<DataType, 3, Frame>,
359 : gr::Tags::SpatialChristoffelSecondKind<DataType, 3, Frame>>;
360 :
361 : template <typename DataType, typename Frame = ::Frame::Inertial>
362 0 : class IntermediateComputer {
363 : public:
364 0 : using CachedBuffer = KerrSchild::CachedBuffer<DataType, Frame>;
365 :
366 0 : IntermediateComputer(const KerrSchild& solution,
367 : const tnsr::I<DataType, 3, Frame>& x);
368 :
369 0 : const KerrSchild& solution() const { return solution_; }
370 :
371 0 : void operator()(
372 : const gsl::not_null<tnsr::I<DataType, 3, Frame>*> x_minus_center,
373 : const gsl::not_null<CachedBuffer*> /*cache*/,
374 : internal_tags::x_minus_center<DataType, Frame> /*meta*/) const;
375 :
376 0 : void operator()(const gsl::not_null<Scalar<DataType>*> a_dot_x,
377 : const gsl::not_null<CachedBuffer*> cache,
378 : internal_tags::a_dot_x<DataType> /*meta*/) const;
379 :
380 0 : void operator()(const gsl::not_null<Scalar<DataType>*> a_dot_x_squared,
381 : const gsl::not_null<CachedBuffer*> cache,
382 : internal_tags::a_dot_x_squared<DataType> /*meta*/) const;
383 :
384 0 : void operator()(const gsl::not_null<Scalar<DataType>*> half_xsq_minus_asq,
385 : const gsl::not_null<CachedBuffer*> cache,
386 : internal_tags::half_xsq_minus_asq<DataType> /*meta*/) const;
387 :
388 0 : void operator()(const gsl::not_null<Scalar<DataType>*> r_squared,
389 : const gsl::not_null<CachedBuffer*> cache,
390 : internal_tags::r_squared<DataType> /*meta*/) const;
391 :
392 0 : void operator()(const gsl::not_null<Scalar<DataType>*> r,
393 : const gsl::not_null<CachedBuffer*> cache,
394 : internal_tags::r<DataType> /*meta*/) const;
395 :
396 0 : void operator()(
397 : const gsl::not_null<Scalar<DataType>*> a_dot_x_over_rsquared,
398 : const gsl::not_null<CachedBuffer*> cache,
399 : internal_tags::a_dot_x_over_rsquared<DataType> /*meta*/) const;
400 :
401 0 : void operator()(const gsl::not_null<Scalar<DataType>*> deriv_log_r_denom,
402 : const gsl::not_null<CachedBuffer*> cache,
403 : internal_tags::deriv_log_r_denom<DataType> /*meta*/) const;
404 :
405 0 : void operator()(
406 : const gsl::not_null<tnsr::i<DataType, 3, Frame>*> deriv_log_r,
407 : const gsl::not_null<CachedBuffer*> cache,
408 : internal_tags::deriv_log_r<DataType, Frame> /*meta*/) const;
409 :
410 0 : void operator()(const gsl::not_null<Scalar<DataType>*> H_denom,
411 : const gsl::not_null<CachedBuffer*> cache,
412 : internal_tags::H_denom<DataType> /*meta*/) const;
413 :
414 0 : void operator()(const gsl::not_null<Scalar<DataType>*> H,
415 : const gsl::not_null<CachedBuffer*> cache,
416 : internal_tags::H<DataType> /*meta*/) const;
417 :
418 0 : void operator()(const gsl::not_null<Scalar<DataType>*> deriv_H_temp1,
419 : const gsl::not_null<CachedBuffer*> cache,
420 : internal_tags::deriv_H_temp1<DataType> /*meta*/) const;
421 :
422 0 : void operator()(const gsl::not_null<Scalar<DataType>*> deriv_H_temp2,
423 : const gsl::not_null<CachedBuffer*> cache,
424 : internal_tags::deriv_H_temp2<DataType> /*meta*/) const;
425 :
426 0 : void operator()(const gsl::not_null<tnsr::i<DataType, 3, Frame>*> deriv_H,
427 : const gsl::not_null<CachedBuffer*> cache,
428 : internal_tags::deriv_H<DataType, Frame> /*meta*/) const;
429 :
430 0 : void operator()(const gsl::not_null<Scalar<DataType>*> denom,
431 : const gsl::not_null<CachedBuffer*> cache,
432 : internal_tags::denom<DataType> /*meta*/) const;
433 :
434 0 : void operator()(const gsl::not_null<Scalar<DataType>*> a_dot_x_over_r,
435 : const gsl::not_null<CachedBuffer*> cache,
436 : internal_tags::a_dot_x_over_r<DataType> /*meta*/) const;
437 :
438 0 : void operator()(const gsl::not_null<tnsr::i<DataType, 3, Frame>*> null_form,
439 : const gsl::not_null<CachedBuffer*> cache,
440 : internal_tags::null_form<DataType, Frame> /*meta*/) const;
441 :
442 0 : void operator()(
443 : const gsl::not_null<tnsr::ij<DataType, 3, Frame>*> deriv_null_form,
444 : const gsl::not_null<CachedBuffer*> cache,
445 : internal_tags::deriv_null_form<DataType, Frame> /*meta*/) const;
446 :
447 0 : void operator()(const gsl::not_null<Scalar<DataType>*> lapse_squared,
448 : const gsl::not_null<CachedBuffer*> cache,
449 : internal_tags::lapse_squared<DataType> /*meta*/) const;
450 :
451 0 : void operator()(const gsl::not_null<Scalar<DataType>*> lapse,
452 : const gsl::not_null<CachedBuffer*> cache,
453 : gr::Tags::Lapse<DataType> /*meta*/) const;
454 :
455 0 : void operator()(
456 : const gsl::not_null<Scalar<DataType>*> deriv_lapse_multiplier,
457 : const gsl::not_null<CachedBuffer*> cache,
458 : internal_tags::deriv_lapse_multiplier<DataType> /*meta*/) const;
459 :
460 0 : void operator()(const gsl::not_null<Scalar<DataType>*> shift_multiplier,
461 : const gsl::not_null<CachedBuffer*> cache,
462 : internal_tags::shift_multiplier<DataType> /*meta*/) const;
463 :
464 0 : void operator()(const gsl::not_null<tnsr::I<DataType, 3, Frame>*> shift,
465 : const gsl::not_null<CachedBuffer*> cache,
466 : gr::Tags::Shift<DataType, 3, Frame> /*meta*/) const;
467 :
468 0 : void operator()(
469 : const gsl::not_null<tnsr::iJ<DataType, 3, Frame>*> deriv_shift,
470 : const gsl::not_null<CachedBuffer*> cache,
471 : DerivShift<DataType, Frame> /*meta*/) const;
472 :
473 0 : void operator()(
474 : const gsl::not_null<tnsr::ii<DataType, 3, Frame>*> spatial_metric,
475 : const gsl::not_null<CachedBuffer*> cache,
476 : gr::Tags::SpatialMetric<DataType, 3, Frame> /*meta*/) const;
477 :
478 0 : void operator()(
479 : const gsl::not_null<tnsr::II<DataType, 3, Frame>*> spatial_metric,
480 : const gsl::not_null<CachedBuffer*> cache,
481 : gr::Tags::InverseSpatialMetric<DataType, 3, Frame> /*meta*/) const;
482 :
483 0 : void operator()(const gsl::not_null<tnsr::ijj<DataType, 3, Frame>*>
484 : deriv_spatial_metric,
485 : const gsl::not_null<CachedBuffer*> cache,
486 : DerivSpatialMetric<DataType, Frame> /*meta*/) const;
487 :
488 0 : void operator()(
489 : const gsl::not_null<tnsr::ii<DataType, 3, Frame>*> dt_spatial_metric,
490 : const gsl::not_null<CachedBuffer*> cache,
491 : ::Tags::dt<gr::Tags::SpatialMetric<DataType, 3, Frame>> /*meta*/) const;
492 :
493 0 : void operator()(
494 : const gsl::not_null<tnsr::ii<DataType, 3, Frame>*> extrinsic_curvature,
495 : const gsl::not_null<CachedBuffer*> cache,
496 : gr::Tags::ExtrinsicCurvature<DataType, 3, Frame> /*meta*/) const;
497 :
498 0 : void operator()(
499 : const gsl::not_null<tnsr::ijj<DataType, 3, Frame>*>
500 : spatial_christoffel_first_kind,
501 : const gsl::not_null<CachedBuffer*> cache,
502 : gr::Tags::SpatialChristoffelFirstKind<DataType, 3, Frame> /*meta*/)
503 : const;
504 :
505 0 : void operator()(
506 : const gsl::not_null<tnsr::Ijj<DataType, 3, Frame>*>
507 : spatial_christoffel_second_kind,
508 : const gsl::not_null<CachedBuffer*> cache,
509 : gr::Tags::SpatialChristoffelSecondKind<DataType, 3, Frame> /*meta*/)
510 : const;
511 :
512 : private:
513 0 : const KerrSchild& solution_;
514 0 : const tnsr::I<DataType, 3, Frame>& x_;
515 : // Here null_vector_0 is simply -1, but if you have a boosted solution,
516 : // then null_vector_0 can be something different, so we leave it coded
517 : // in instead of eliminating it.
518 0 : static constexpr double null_vector_0_ = -1.0;
519 : };
520 :
521 : template <typename DataType, typename Frame = ::Frame::Inertial>
522 0 : class IntermediateVars : public CachedBuffer<DataType, Frame> {
523 : public:
524 0 : using CachedBuffer = KerrSchild::CachedBuffer<DataType, Frame>;
525 : using CachedBuffer::CachedBuffer;
526 1 : using CachedBuffer::get_var;
527 :
528 0 : tnsr::i<DataType, 3, Frame> get_var(
529 : const IntermediateComputer<DataType, Frame>& computer,
530 : DerivLapse<DataType, Frame> /*meta*/);
531 :
532 0 : Scalar<DataType> get_var(
533 : const IntermediateComputer<DataType, Frame>& computer,
534 : ::Tags::dt<gr::Tags::Lapse<DataType>> /*meta*/);
535 :
536 0 : tnsr::I<DataType, 3, Frame> get_var(
537 : const IntermediateComputer<DataType, Frame>& computer,
538 : ::Tags::dt<gr::Tags::Shift<DataType, 3, Frame>> /*meta*/);
539 :
540 0 : Scalar<DataType> get_var(
541 : const IntermediateComputer<DataType, Frame>& computer,
542 : gr::Tags::SqrtDetSpatialMetric<DataType> /*meta*/);
543 :
544 0 : tnsr::i<DataType, 3, Frame> get_var(
545 : const IntermediateComputer<DataType, Frame>& computer,
546 : gr::Tags::DerivDetSpatialMetric<DataType, 3, Frame> /*meta*/);
547 :
548 0 : Scalar<DataType> get_var(
549 : const IntermediateComputer<DataType, Frame>& computer,
550 : gr::Tags::TraceExtrinsicCurvature<DataType> /*meta*/);
551 :
552 0 : tnsr::I<DataType, 3, Frame> get_var(
553 : const IntermediateComputer<DataType, Frame>& computer,
554 : gr::Tags::TraceSpatialChristoffelSecondKind<DataType, 3,
555 : Frame> /*meta*/);
556 0 : tnsr::Abb<DataType, 3, Frame> get_var(
557 : const IntermediateComputer<DataType, Frame>& computer,
558 : gr::Tags::SpacetimeChristoffelSecondKind<DataType, 3, Frame> /*meta*/);
559 :
560 : private:
561 : // Here null_vector_0 is simply -1, but if you have a boosted solution,
562 : // then null_vector_0 can be something different, so we leave it coded
563 : // in instead of eliminating it.
564 0 : static constexpr double null_vector_0_ = -1.0;
565 : };
566 :
567 : private:
568 0 : double mass_{std::numeric_limits<double>::signaling_NaN()};
569 0 : std::array<double, volume_dim> dimensionless_spin_ =
570 : make_array<volume_dim>(std::numeric_limits<double>::signaling_NaN());
571 0 : std::array<double, volume_dim> center_ =
572 : make_array<volume_dim>(std::numeric_limits<double>::signaling_NaN());
573 0 : bool zero_spin_{};
574 : };
575 :
576 0 : SPECTRE_ALWAYS_INLINE bool operator==(const KerrSchild& lhs,
577 : const KerrSchild& rhs) {
578 : return lhs.mass() == rhs.mass() and
579 : lhs.dimensionless_spin() == rhs.dimensionless_spin() and
580 : lhs.center() == rhs.center();
581 : }
582 :
583 0 : SPECTRE_ALWAYS_INLINE bool operator!=(const KerrSchild& lhs,
584 : const KerrSchild& rhs) {
585 : return not(lhs == rhs);
586 : }
587 : } // namespace Solutions
588 : } // namespace gr
|
