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/TaggedTuple.hpp"
12 : #include "DataStructures/Tags/TempTensor.hpp"
13 : #include "DataStructures/Tensor/TypeAliases.hpp"
14 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
15 : #include "Options/Context.hpp"
16 : #include "Options/String.hpp"
17 : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
18 : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Solutions.hpp"
19 : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"
20 : #include "Utilities/ContainerHelpers.hpp"
21 : #include "Utilities/ForceInline.hpp"
22 : #include "Utilities/MakeArray.hpp"
23 : #include "Utilities/TMPL.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 : * \gamma_{i j} &=& \delta_{i j} + 2 H l_i l_j,\\
81 : * \gamma^{i j} &=& \delta^{i j} - {2 H l^i l^j \over 1+2H l^t l^t},\\
82 : * {\rm det} \gamma_{i j}&=& 1+2H l^t l^t,\\
83 : * \partial_k ({\rm det} \gamma_{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}\gamma_{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$\gamma^{i j} l_j\f$; it is equal
116 : * to \f$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 : * ## Boost of the Kerr-Schild solution
219 : *
220 : * We add initial momentum to the solution by applying a Lorentz boost to the
221 : * metric. Since the Kerr-Schild metric can be expressed covariantly in terms of
222 : * the Minkowski metric, a scalar function and a one form, we construct the
223 : * metric in the rest frame of the black hole and then apply an inverse boost to
224 : * each of the covariant objects individually. Notice that we also need to
225 : * appropriately boost the coordinates to the to the rest frame before computing
226 : * the metric.
227 : *
228 : * \warning While technically the boosted Kerr-Schild metric is dependent on
229 : * both the time and space coordinates, we have implemented it only at $t = 0$
230 : * as in SpEC. Therefore it is technically not an analytic solution and should
231 : * not be used to compute errors with respect to it.
232 : *
233 : * Moreover, since the boosted solution is intended for use as initial data,
234 : * we do not compute the time derivatives of the lapse and shift in the boosted
235 : * frame but set them to zero.
236 : *
237 : * Consequently, the gr::Tags::SpacetimeChristoffelSecondKind computed here,
238 : * corresponds to the boosted Kerr-Schild for the gauge where lapse and shift
239 : * have vanishing derivatives.
240 : *
241 : */
242 1 : class KerrSchild : public AnalyticSolution<3_st>,
243 : public MarkAsAnalyticSolution {
244 : public:
245 0 : struct Mass {
246 0 : using type = double;
247 0 : static constexpr Options::String help = {"Mass of the black hole"};
248 0 : static type lower_bound() { return 0.; }
249 : };
250 0 : struct Spin {
251 0 : using type = std::array<double, volume_dim>;
252 0 : static constexpr Options::String help = {
253 : "The [x,y,z] dimensionless spin of the black hole"};
254 : };
255 0 : struct Center {
256 0 : using type = std::array<double, volume_dim>;
257 0 : static constexpr Options::String help = {
258 : "The [x,y,z] center of the black hole"};
259 : };
260 0 : struct Velocity {
261 0 : using type = std::array<double, volume_dim>;
262 0 : static constexpr Options::String help = {
263 : "The [x,y,z] boost velocity of the black hole"};
264 : };
265 0 : using options = tmpl::list<Mass, Spin, Center, Velocity>;
266 0 : static constexpr Options::String help{
267 : "Black hole in Kerr-Schild coordinates"};
268 :
269 0 : KerrSchild(double mass, const std::array<double, 3>& dimensionless_spin,
270 : const std::array<double, 3>& center,
271 : const std::array<double, 3>& boost_velocity = {{0., 0., 0.}},
272 : const Options::Context& context = {});
273 :
274 0 : explicit KerrSchild(CkMigrateMessage* /*msg*/);
275 :
276 : template <typename DataType, typename Frame = Frame::Inertial>
277 0 : using tags = tmpl::flatten<tmpl::list<
278 : AnalyticSolution<3_st>::tags<DataType, Frame>,
279 : gr::Tags::DerivDetSpatialMetric<DataType, 3, Frame>,
280 : gr::Tags::TraceExtrinsicCurvature<DataType>,
281 : gr::Tags::SpatialChristoffelFirstKind<DataType, 3, Frame>,
282 : gr::Tags::SpatialChristoffelSecondKind<DataType, 3, Frame>,
283 : gr::Tags::TraceSpatialChristoffelSecondKind<DataType, 3, Frame>,
284 : gr::Tags::SpacetimeChristoffelSecondKind<DataType, 3, Frame>>>;
285 :
286 0 : KerrSchild() = default;
287 0 : KerrSchild(const KerrSchild& /*rhs*/) = default;
288 0 : KerrSchild& operator=(const KerrSchild& /*rhs*/) = default;
289 0 : KerrSchild(KerrSchild&& /*rhs*/) = default;
290 0 : KerrSchild& operator=(KerrSchild&& /*rhs*/) = default;
291 0 : ~KerrSchild() = default;
292 :
293 : template <typename DataType, typename Frame, typename... Tags>
294 0 : tuples::TaggedTuple<Tags...> variables(
295 : const tnsr::I<DataType, volume_dim, Frame>& x, double /*t*/,
296 : tmpl::list<Tags...> /*meta*/) const {
297 : static_assert(
298 : tmpl2::flat_all_v<
299 : tmpl::list_contains_v<tags<DataType, Frame>, Tags>...>,
300 : "At least one of the requested tags is not supported. The requested "
301 : "tags are listed as template parameters of the `variables` function.");
302 : IntermediateVars<DataType, Frame> cache(get_size(*x.begin()));
303 : IntermediateComputer<DataType, Frame> computer(*this, x);
304 : return {cache.get_var(computer, Tags{})...};
305 : }
306 :
307 : // NOLINTNEXTLINE(google-runtime-references)
308 0 : void pup(PUP::er& p);
309 :
310 0 : double mass() const { return mass_; }
311 0 : const std::array<double, volume_dim>& center() const { return center_; }
312 0 : const std::array<double, volume_dim>& dimensionless_spin() const {
313 : return dimensionless_spin_;
314 : }
315 0 : const std::array<double, volume_dim>& boost_velocity() const {
316 : return boost_velocity_;
317 : }
318 0 : bool zero_spin() const { return zero_spin_; }
319 0 : bool zero_velocity() const { return zero_velocity_; }
320 :
321 0 : struct internal_tags {
322 : template <typename DataType, typename Frame = ::Frame::Inertial>
323 0 : using x_minus_center_unboosted = ::Tags::TempI<0, 3, Frame, DataType>;
324 : template <typename DataType, typename Frame = ::Frame::Inertial>
325 0 : using x_minus_center = ::Tags::TempI<1, 3, Frame, DataType>;
326 : template <typename DataType>
327 0 : using a_dot_x = ::Tags::TempScalar<2, DataType>;
328 : template <typename DataType>
329 0 : using a_dot_x_squared = ::Tags::TempScalar<3, DataType>;
330 : template <typename DataType>
331 0 : using half_xsq_minus_asq = ::Tags::TempScalar<4, DataType>;
332 : template <typename DataType>
333 0 : using r_squared = ::Tags::TempScalar<5, DataType>;
334 : template <typename DataType>
335 0 : using r = ::Tags::TempScalar<6, DataType>;
336 : template <typename DataType>
337 0 : using a_dot_x_over_rsquared = ::Tags::TempScalar<7, DataType>;
338 : template <typename DataType>
339 0 : using deriv_log_r_denom = ::Tags::TempScalar<8, DataType>;
340 : template <typename DataType, typename Frame = ::Frame::Inertial>
341 0 : using deriv_log_r = ::Tags::Tempi<9, 3, Frame, DataType>;
342 : template <typename DataType>
343 0 : using H_denom = ::Tags::TempScalar<10, DataType>;
344 : template <typename DataType>
345 0 : using H = ::Tags::TempScalar<11, DataType>;
346 : template <typename DataType>
347 0 : using deriv_H_temp1 = ::Tags::TempScalar<12, DataType>;
348 : template <typename DataType>
349 0 : using deriv_H_temp2 = ::Tags::TempScalar<13, DataType>;
350 : template <typename DataType, typename Frame = ::Frame::Inertial>
351 0 : using deriv_H_unboosted = ::Tags::Tempa<14, 3, Frame, DataType>;
352 : template <typename DataType, typename Frame = ::Frame::Inertial>
353 0 : using deriv_H = ::Tags::Tempa<15, 3, Frame, DataType>;
354 : template <typename DataType>
355 0 : using denom = ::Tags::TempScalar<16, DataType>;
356 : template <typename DataType>
357 0 : using a_dot_x_over_r = ::Tags::TempScalar<17, DataType>;
358 : template <typename DataType, typename Frame = ::Frame::Inertial>
359 0 : using null_form_unboosted = ::Tags::Tempa<18, 3, Frame, DataType>;
360 : template <typename DataType, typename Frame = ::Frame::Inertial>
361 0 : using null_form = ::Tags::Tempa<19, 3, Frame, DataType>;
362 : template <typename DataType, typename Frame = ::Frame::Inertial>
363 0 : using deriv_null_form_unboosted = ::Tags::Tempab<20, 3, Frame, DataType>;
364 : template <typename DataType, typename Frame = ::Frame::Inertial>
365 0 : using deriv_null_form = ::Tags::Tempab<21, 3, Frame, DataType>;
366 : template <typename DataType>
367 0 : using null_form_dot_deriv_H = ::Tags::TempScalar<22, DataType>;
368 : template <typename DataType, typename Frame = ::Frame::Inertial>
369 0 : using null_form_dot_deriv_null_form = ::Tags::Tempi<23, 3, Frame, DataType>;
370 : template <typename DataType>
371 0 : using lapse_squared = ::Tags::TempScalar<24, DataType>;
372 : template <typename DataType>
373 0 : using deriv_lapse_multiplier = ::Tags::TempScalar<25, DataType>;
374 : template <typename DataType>
375 0 : using shift_multiplier = ::Tags::TempScalar<26, DataType>;
376 : };
377 :
378 : template <typename DataType, typename Frame = ::Frame::Inertial>
379 0 : using CachedBuffer = CachedTempBuffer<
380 : internal_tags::x_minus_center_unboosted<DataType, Frame>,
381 : internal_tags::x_minus_center<DataType, Frame>,
382 : internal_tags::a_dot_x<DataType>,
383 : internal_tags::a_dot_x_squared<DataType>,
384 : internal_tags::half_xsq_minus_asq<DataType>,
385 : internal_tags::r_squared<DataType>, internal_tags::r<DataType>,
386 : internal_tags::a_dot_x_over_rsquared<DataType>,
387 : internal_tags::deriv_log_r_denom<DataType>,
388 : internal_tags::deriv_log_r<DataType, Frame>,
389 : internal_tags::H_denom<DataType>, internal_tags::H<DataType>,
390 : internal_tags::deriv_H_temp1<DataType>,
391 : internal_tags::deriv_H_temp2<DataType>,
392 : internal_tags::deriv_H_unboosted<DataType, Frame>,
393 : internal_tags::deriv_H<DataType, Frame>, internal_tags::denom<DataType>,
394 : internal_tags::a_dot_x_over_r<DataType>,
395 : internal_tags::null_form_unboosted<DataType, Frame>,
396 : internal_tags::null_form<DataType, Frame>,
397 : internal_tags::deriv_null_form_unboosted<DataType, Frame>,
398 : internal_tags::deriv_null_form<DataType, Frame>,
399 : internal_tags::null_form_dot_deriv_H<DataType>,
400 : internal_tags::null_form_dot_deriv_null_form<DataType, Frame>,
401 : internal_tags::lapse_squared<DataType>, gr::Tags::Lapse<DataType>,
402 : internal_tags::deriv_lapse_multiplier<DataType>,
403 : internal_tags::shift_multiplier<DataType>,
404 : gr::Tags::Shift<DataType, 3, Frame>, DerivShift<DataType, Frame>,
405 : gr::Tags::SpatialMetric<DataType, 3, Frame>,
406 : gr::Tags::InverseSpatialMetric<DataType, 3, Frame>,
407 : DerivSpatialMetric<DataType, Frame>,
408 : ::Tags::dt<gr::Tags::SpatialMetric<DataType, 3, Frame>>,
409 : gr::Tags::ExtrinsicCurvature<DataType, 3, Frame>,
410 : gr::Tags::SpatialChristoffelFirstKind<DataType, 3, Frame>,
411 : gr::Tags::SpatialChristoffelSecondKind<DataType, 3, Frame>>;
412 :
413 : template <typename DataType, typename Frame = ::Frame::Inertial>
414 0 : class IntermediateComputer {
415 : public:
416 0 : using CachedBuffer = KerrSchild::CachedBuffer<DataType, Frame>;
417 :
418 0 : IntermediateComputer(const KerrSchild& solution,
419 : const tnsr::I<DataType, 3, Frame>& x);
420 :
421 0 : const KerrSchild& solution() const { return solution_; }
422 :
423 0 : void operator()(
424 : gsl::not_null<tnsr::I<DataType, 3, Frame>*> x_minus_center,
425 : gsl::not_null<CachedBuffer*> /*cache*/,
426 : internal_tags::x_minus_center_unboosted<DataType, Frame> /*meta*/)
427 : const;
428 :
429 0 : void operator()(
430 : gsl::not_null<tnsr::I<DataType, 3, Frame>*> x_minus_center_boosted,
431 : gsl::not_null<CachedBuffer*> /*cache*/,
432 : internal_tags::x_minus_center<DataType, Frame> /*meta*/) const;
433 :
434 0 : void operator()(gsl::not_null<Scalar<DataType>*> a_dot_x,
435 : gsl::not_null<CachedBuffer*> cache,
436 : internal_tags::a_dot_x<DataType> /*meta*/) const;
437 :
438 0 : void operator()(gsl::not_null<Scalar<DataType>*> a_dot_x_squared,
439 : gsl::not_null<CachedBuffer*> cache,
440 : internal_tags::a_dot_x_squared<DataType> /*meta*/) const;
441 :
442 0 : void operator()(gsl::not_null<Scalar<DataType>*> half_xsq_minus_asq,
443 : gsl::not_null<CachedBuffer*> cache,
444 : internal_tags::half_xsq_minus_asq<DataType> /*meta*/) const;
445 :
446 0 : void operator()(gsl::not_null<Scalar<DataType>*> r_squared,
447 : gsl::not_null<CachedBuffer*> cache,
448 : internal_tags::r_squared<DataType> /*meta*/) const;
449 :
450 0 : void operator()(gsl::not_null<Scalar<DataType>*> r,
451 : gsl::not_null<CachedBuffer*> cache,
452 : internal_tags::r<DataType> /*meta*/) const;
453 :
454 0 : void operator()(
455 : gsl::not_null<Scalar<DataType>*> a_dot_x_over_rsquared,
456 : gsl::not_null<CachedBuffer*> cache,
457 : internal_tags::a_dot_x_over_rsquared<DataType> /*meta*/) const;
458 :
459 0 : void operator()(gsl::not_null<Scalar<DataType>*> deriv_log_r_denom,
460 : gsl::not_null<CachedBuffer*> cache,
461 : internal_tags::deriv_log_r_denom<DataType> /*meta*/) const;
462 :
463 0 : void operator()(gsl::not_null<tnsr::i<DataType, 3, Frame>*> deriv_log_r,
464 : gsl::not_null<CachedBuffer*> cache,
465 : internal_tags::deriv_log_r<DataType, Frame> /*meta*/) const;
466 :
467 0 : void operator()(gsl::not_null<Scalar<DataType>*> H_denom,
468 : gsl::not_null<CachedBuffer*> cache,
469 : internal_tags::H_denom<DataType> /*meta*/) const;
470 :
471 0 : void operator()(gsl::not_null<Scalar<DataType>*> H,
472 : gsl::not_null<CachedBuffer*> cache,
473 : internal_tags::H<DataType> /*meta*/) const;
474 :
475 0 : void operator()(gsl::not_null<Scalar<DataType>*> deriv_H_temp1,
476 : gsl::not_null<CachedBuffer*> cache,
477 : internal_tags::deriv_H_temp1<DataType> /*meta*/) const;
478 :
479 0 : void operator()(gsl::not_null<Scalar<DataType>*> deriv_H_temp2,
480 : gsl::not_null<CachedBuffer*> cache,
481 : internal_tags::deriv_H_temp2<DataType> /*meta*/) const;
482 :
483 0 : void operator()(
484 : gsl::not_null<tnsr::a<DataType, 3, Frame>*> deriv_H,
485 : gsl::not_null<CachedBuffer*> cache,
486 : internal_tags::deriv_H_unboosted<DataType, Frame> /*meta*/) const;
487 :
488 0 : void operator()(gsl::not_null<tnsr::a<DataType, 3, Frame>*> deriv_H_boosted,
489 : gsl::not_null<CachedBuffer*> cache,
490 : internal_tags::deriv_H<DataType, Frame> /*meta*/) const;
491 :
492 0 : void operator()(gsl::not_null<Scalar<DataType>*> denom,
493 : gsl::not_null<CachedBuffer*> cache,
494 : internal_tags::denom<DataType> /*meta*/) const;
495 :
496 0 : void operator()(gsl::not_null<Scalar<DataType>*> a_dot_x_over_r,
497 : gsl::not_null<CachedBuffer*> cache,
498 : internal_tags::a_dot_x_over_r<DataType> /*meta*/) const;
499 :
500 0 : void operator()(
501 : gsl::not_null<tnsr::a<DataType, 3, Frame>*> null_form,
502 : gsl::not_null<CachedBuffer*> cache,
503 : internal_tags::null_form_unboosted<DataType, Frame> /*meta*/) const;
504 :
505 0 : void operator()(
506 : gsl::not_null<tnsr::a<DataType, 3, Frame>*> null_form_boosted,
507 : gsl::not_null<CachedBuffer*> cache,
508 : internal_tags::null_form<DataType, Frame> /*meta*/) const;
509 :
510 0 : void operator()(
511 : gsl::not_null<tnsr::ab<DataType, 3, Frame>*> deriv_null_form,
512 : gsl::not_null<CachedBuffer*> cache,
513 : internal_tags::deriv_null_form_unboosted<DataType, Frame> /*meta*/)
514 : const;
515 :
516 0 : void operator()(
517 : gsl::not_null<tnsr::ab<DataType, 3, Frame>*> deriv_null_form_boosted,
518 : gsl::not_null<CachedBuffer*> cache,
519 : internal_tags::deriv_null_form<DataType, Frame> /*meta*/) const;
520 :
521 0 : void operator()(gsl::not_null<Scalar<DataType>*> lapse_squared,
522 : gsl::not_null<CachedBuffer*> cache,
523 : internal_tags::lapse_squared<DataType> /*meta*/) const;
524 :
525 0 : void operator()(gsl::not_null<Scalar<DataType>*> lapse,
526 : gsl::not_null<CachedBuffer*> cache,
527 : gr::Tags::Lapse<DataType> /*meta*/) const;
528 :
529 0 : void operator()(
530 : gsl::not_null<Scalar<DataType>*> deriv_lapse_multiplier,
531 : gsl::not_null<CachedBuffer*> cache,
532 : internal_tags::deriv_lapse_multiplier<DataType> /*meta*/) const;
533 :
534 0 : void operator()(gsl::not_null<Scalar<DataType>*> shift_multiplier,
535 : gsl::not_null<CachedBuffer*> cache,
536 : internal_tags::shift_multiplier<DataType> /*meta*/) const;
537 :
538 0 : void operator()(gsl::not_null<tnsr::I<DataType, 3, Frame>*> shift,
539 : gsl::not_null<CachedBuffer*> cache,
540 : gr::Tags::Shift<DataType, 3, Frame> /*meta*/) const;
541 :
542 0 : void operator()(gsl::not_null<tnsr::iJ<DataType, 3, Frame>*> deriv_shift,
543 : gsl::not_null<CachedBuffer*> cache,
544 : DerivShift<DataType, Frame> /*meta*/) const;
545 :
546 0 : void operator()(gsl::not_null<tnsr::ii<DataType, 3, Frame>*> spatial_metric,
547 : gsl::not_null<CachedBuffer*> cache,
548 : gr::Tags::SpatialMetric<DataType, 3, Frame> /*meta*/) const;
549 :
550 0 : void operator()(
551 : gsl::not_null<tnsr::II<DataType, 3, Frame>*> spatial_metric,
552 : gsl::not_null<CachedBuffer*> cache,
553 : gr::Tags::InverseSpatialMetric<DataType, 3, Frame> /*meta*/) const;
554 :
555 0 : void operator()(
556 : gsl::not_null<tnsr::ijj<DataType, 3, Frame>*> deriv_spatial_metric,
557 : gsl::not_null<CachedBuffer*> cache,
558 : DerivSpatialMetric<DataType, Frame> /*meta*/) const;
559 :
560 0 : void operator()(
561 : gsl::not_null<tnsr::ii<DataType, 3, Frame>*> dt_spatial_metric,
562 : gsl::not_null<CachedBuffer*> cache,
563 : ::Tags::dt<gr::Tags::SpatialMetric<DataType, 3, Frame>> /*meta*/) const;
564 :
565 0 : void operator()(
566 : gsl::not_null<Scalar<DataType>*> null_form_dot_deriv_H,
567 : gsl::not_null<CachedBuffer*> cache,
568 : internal_tags::null_form_dot_deriv_H<DataType> /*meta*/) const;
569 :
570 0 : void operator()(
571 : gsl::not_null<tnsr::i<DataType, 3, Frame>*>
572 : null_form_dot_deriv_null_form,
573 : gsl::not_null<CachedBuffer*> cache,
574 : internal_tags::null_form_dot_deriv_null_form<DataType, Frame> /*meta*/)
575 : const;
576 :
577 0 : void operator()(
578 : gsl::not_null<tnsr::ii<DataType, 3, Frame>*> extrinsic_curvature,
579 : gsl::not_null<CachedBuffer*> cache,
580 : gr::Tags::ExtrinsicCurvature<DataType, 3, Frame> /*meta*/) const;
581 :
582 0 : void operator()(
583 : gsl::not_null<tnsr::ijj<DataType, 3, Frame>*>
584 : spatial_christoffel_first_kind,
585 : gsl::not_null<CachedBuffer*> cache,
586 : gr::Tags::SpatialChristoffelFirstKind<DataType, 3, Frame> /*meta*/)
587 : const;
588 :
589 0 : void operator()(
590 : gsl::not_null<tnsr::Ijj<DataType, 3, Frame>*>
591 : spatial_christoffel_second_kind,
592 : gsl::not_null<CachedBuffer*> cache,
593 : gr::Tags::SpatialChristoffelSecondKind<DataType, 3, Frame> /*meta*/)
594 : const;
595 :
596 : private:
597 0 : const KerrSchild& solution_;
598 0 : const tnsr::I<DataType, 3, Frame>& x_;
599 : // Here null_vector_0 is simply -1, but if you have a boosted solution,
600 : // then null_vector_0 can be something different, so we leave it coded
601 : // in instead of eliminating it.
602 0 : static constexpr double null_vector_0_ = -1.0;
603 : };
604 :
605 : template <typename DataType, typename Frame = ::Frame::Inertial>
606 0 : class IntermediateVars : public CachedBuffer<DataType, Frame> {
607 : public:
608 0 : using CachedBuffer = KerrSchild::CachedBuffer<DataType, Frame>;
609 : using CachedBuffer::CachedBuffer;
610 1 : using CachedBuffer::get_var;
611 :
612 0 : tnsr::i<DataType, 3, Frame> get_var(
613 : const IntermediateComputer<DataType, Frame>& computer,
614 : DerivLapse<DataType, Frame> /*meta*/);
615 :
616 0 : Scalar<DataType> get_var(
617 : const IntermediateComputer<DataType, Frame>& computer,
618 : ::Tags::dt<gr::Tags::Lapse<DataType>> /*meta*/);
619 :
620 0 : tnsr::I<DataType, 3, Frame> get_var(
621 : const IntermediateComputer<DataType, Frame>& computer,
622 : ::Tags::dt<gr::Tags::Shift<DataType, 3, Frame>> /*meta*/);
623 :
624 0 : Scalar<DataType> get_var(
625 : const IntermediateComputer<DataType, Frame>& computer,
626 : gr::Tags::SqrtDetSpatialMetric<DataType> /*meta*/);
627 :
628 0 : tnsr::i<DataType, 3, Frame> get_var(
629 : const IntermediateComputer<DataType, Frame>& computer,
630 : gr::Tags::DerivDetSpatialMetric<DataType, 3, Frame> /*meta*/);
631 :
632 0 : Scalar<DataType> get_var(
633 : const IntermediateComputer<DataType, Frame>& computer,
634 : gr::Tags::TraceExtrinsicCurvature<DataType> /*meta*/);
635 :
636 0 : tnsr::I<DataType, 3, Frame> get_var(
637 : const IntermediateComputer<DataType, Frame>& computer,
638 : gr::Tags::TraceSpatialChristoffelSecondKind<DataType, 3,
639 : Frame> /*meta*/);
640 0 : tnsr::Abb<DataType, 3, Frame> get_var(
641 : const IntermediateComputer<DataType, Frame>& computer,
642 : gr::Tags::SpacetimeChristoffelSecondKind<DataType, 3, Frame> /*meta*/);
643 :
644 : private:
645 : // Here null_vector_0 is simply -1, but if you have a boosted solution,
646 : // then null_vector_0 can be something different, so we leave it coded
647 : // in instead of eliminating it.
648 0 : static constexpr double null_vector_0_ = -1.0;
649 : };
650 :
651 : private:
652 0 : double mass_{std::numeric_limits<double>::signaling_NaN()};
653 0 : std::array<double, volume_dim> dimensionless_spin_ =
654 : make_array<volume_dim>(std::numeric_limits<double>::signaling_NaN());
655 0 : std::array<double, volume_dim> center_ =
656 : make_array<volume_dim>(std::numeric_limits<double>::signaling_NaN());
657 0 : std::array<double, volume_dim> boost_velocity_ =
658 : make_array<volume_dim>(std::numeric_limits<double>::signaling_NaN());
659 0 : bool zero_spin_{};
660 0 : bool zero_velocity_{};
661 : };
662 :
663 0 : SPECTRE_ALWAYS_INLINE bool operator==(const KerrSchild& lhs,
664 : const KerrSchild& rhs) {
665 : return lhs.mass() == rhs.mass() and
666 : lhs.dimensionless_spin() == rhs.dimensionless_spin() and
667 : lhs.center() == rhs.center() and
668 : lhs.boost_velocity() == rhs.boost_velocity();
669 : }
670 :
671 0 : SPECTRE_ALWAYS_INLINE bool operator!=(const KerrSchild& lhs,
672 : const KerrSchild& rhs) {
673 : return not(lhs == rhs);
674 : }
675 : } // namespace Solutions
676 : } // namespace gr
|
