Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <complex>
7 : #include <cstddef>
8 : #include <limits>
9 : #include <memory>
10 :
11 : #include "DataStructures/DataBox/Prefixes.hpp"
12 : #include "DataStructures/SpinWeighted.hpp"
13 : #include "DataStructures/Tensor/Tensor.hpp"
14 : #include "Evolution/Systems/Cce/AnalyticSolutions/SphericalMetricData.hpp"
15 : #include "Evolution/Systems/Cce/AnalyticSolutions/WorldtubeData.hpp"
16 : #include "Evolution/Systems/Cce/Tags.hpp"
17 : #include "Options/String.hpp"
18 : #include "Utilities/Gsl.hpp"
19 : #include "Utilities/Serialization/CharmPupable.hpp"
20 : #include "Utilities/TMPL.hpp"
21 :
22 : /// \cond
23 : class DataVector;
24 : class ComplexDataVector;
25 : /// \endcond
26 :
27 : namespace Cce {
28 : namespace Solutions {
29 :
30 : /*!
31 : * \brief Computes the analytic data for a Teukolsky wave solution described in
32 : * \cite Barkett2019uae.
33 : *
34 : * \details This test computes an outgoing perturbative wave solution in
35 : * spherical coordinates with wave profile
36 : *
37 : * \f[
38 : * F(u) = A e^{- u^2 / k^2}.
39 : * \f]
40 : */
41 1 : struct TeukolskyWave : public SphericalMetricData {
42 0 : struct ExtractionRadius {
43 0 : using type = double;
44 0 : static constexpr Options::String help{
45 : "The extraction radius of the spherical solution"};
46 0 : static type lower_bound() { return 0.0; }
47 : };
48 0 : struct Amplitude {
49 0 : using type = double;
50 0 : static constexpr Options::String help{
51 : "The amplitude of the Teukolsky wave."};
52 0 : static type lower_bound() { return 0.0; }
53 : };
54 0 : struct Duration {
55 0 : using type = double;
56 0 : static constexpr Options::String help{
57 : "The characteristic duration of the Gaussian envelope."};
58 0 : static type lower_bound() { return 0.0; }
59 : };
60 :
61 0 : static constexpr Options::String help{
62 : "An analytic solution derived from the linearized Teukolsky equation"};
63 :
64 0 : using options = tmpl::list<ExtractionRadius, Amplitude, Duration>;
65 :
66 0 : WRAPPED_PUPable_decl_template(TeukolskyWave); // NOLINT
67 :
68 0 : explicit TeukolskyWave(CkMigrateMessage* msg) : SphericalMetricData(msg) {}
69 :
70 : // clang doesn't manage to use = default correctly in this case
71 : // NOLINTNEXTLINE(modernize-use-equals-default)
72 0 : TeukolskyWave() {}
73 :
74 0 : TeukolskyWave(double extraction_radius, double amplitude, double duration);
75 :
76 0 : std::unique_ptr<WorldtubeData> get_clone() const override;
77 :
78 0 : void pup(PUP::er& p) override;
79 :
80 : private:
81 : /*
82 : * The A coefficient is calculated as
83 : *
84 : * \f[
85 : * A = 3 (3 k^4 + 4 r^2 u^2 - 2 k^2 r^2 (r + 3 u)) a e^{-u^2/k^2} / (k^4 r^5)
86 : * \f]
87 : *
88 : * where \f$r\f$ is the extraction radius, \f$u\f$ the retarded time, \f$a\f$
89 : * the amplitude of the pulse and \f$k\f$ the duration of the pulse.
90 : */
91 0 : double pulse_profile_coefficient_a(double time) const;
92 :
93 : /*
94 : * The B coefficient is calculated as
95 : *
96 : * \f[
97 : * B = 2 (-3 k^6 + 4 r^4 u^3 - 6 k^2 r^3 u (r + u) + 3 k^4 r^2 (r + 2 u))
98 : * * a e^{-u^2/k^2} / (k^6 r^6)
99 : * \f]
100 : *
101 : * where \f$r\f$ is the extraction radius, \f$u\f$ the retarded time, \f$a\f$
102 : * the amplitude of the pulse and \f$k\f$ the duration of the pulse.
103 : */
104 0 : double pulse_profile_coefficient_b(double time) const;
105 :
106 : /*
107 : * The C coefficient is calculated as
108 : *
109 : * \f[
110 : * C = (1/4) (21 k^8 + 16 r^4 u^4 - 16 k^2 r^3 u^2 (3 r + u)
111 : * - 6 k^6 r (3 r + 7 u) + 12 k^4 r^2 (r^2 + 2 r u + 3 u^2))
112 : * * a e^{-u^2/k^2} / (k^8 r^5)
113 : * \f]
114 : *
115 : * where \f$r\f$ is the extraction radius, \f$u\f$ the retarded time, \f$a\f$
116 : * the amplitude of the pulse and \f$k\f$ the duration of the pulse.
117 : */
118 0 : double pulse_profile_coefficient_c(double time) const;
119 :
120 : /*
121 : * The time derivative of the A coefficient is calculated as
122 : *
123 : * \f[
124 : * \partial_t A = -6 (4 r^2 u^3 + 3 k^4 (r + u) - 6 k^2 r u (r + u))
125 : * a e^{-u^2/k^2} / (k^6 r^5)
126 : * \f]
127 : *
128 : * where \f$r\f$ is the extraction radius, \f$u\f$ the retarded time, \f$a\f$
129 : * the amplitude of the pulse and \f$k\f$ the duration of the pulse.
130 : */
131 0 : double dt_pulse_profile_coefficient_a(double time) const;
132 :
133 : /*
134 : * The time derivative of the B coefficient is calculated as
135 : *
136 : * \f[
137 : * \partial_t B = 4 (-4 r^4 u^4 + 6 k ^2 r^3 u^2 (2 r + u)
138 : * + 3 k^6 (r^2 + u) - 3 k^4 r^2 (r + u) (r + 2u))
139 : * a e^{-u^2/k^2} / (k^8 r^6)
140 : * \f]
141 : *
142 : * where \f$r\f$ is the extraction radius, \f$u\f$ the retarded time, \f$a\f$
143 : * the amplitude of the pulse and \f$k\f$ the duration of the pulse.
144 : */
145 0 : double dt_pulse_profile_coefficient_b(double time) const;
146 :
147 : /*
148 : * The time derivative of the C coefficient is calculated as
149 : *
150 : * \f[
151 : * \partial_t C = (1/2) (-16 r^4 u^5 - 21 k^8 (r + u)
152 : * + 16 k^2 r^3 u^3 (5 r + u) + 6 k^6 r (r + u) (2 r + 7 u)
153 : * - 12 k^4 r^2 u (4 r^2 + 4 r u + 3 u^2))
154 : * a e^{-u^2/k^2} / (k^10 r^5)
155 : * \f]
156 : *
157 : * where \f$r\f$ is the extraction radius, \f$u\f$ the retarded time, \f$a\f$
158 : * the amplitude of the pulse and \f$k\f$ the duration of the pulse.
159 : */
160 0 : double dt_pulse_profile_coefficient_c(double time) const;
161 :
162 : /*
163 : * The radial derivative of the A coefficient is calculated as
164 : *
165 : * \f[
166 : * \partial_r A + \partial_t A = - 9 (5 k^4 + 4 r^2 u^2 - 2 k^2 r (r + 4 u))
167 : * a e^{-u^2/k^2} / (k^4 r^6)
168 : * \f]
169 : *
170 : * where \f$r\f$ is the extraction radius, \f$u\f$ the retarded time, \f$a\f$
171 : * the amplitude of the pulse and \f$k\f$ the duration of the pulse.
172 : */
173 0 : double dr_pulse_profile_coefficient_a(double time) const;
174 :
175 : /*
176 : * The radial derivative of the B coefficient is calculated as
177 : *
178 : * \f[
179 : * \partial_r B + \partial_t B = 2 (18 k^6 - 8 r^4 u^3 + 6 k^2 r^3 u (2 r + 3
180 : * u)
181 : * - 3 k^4 r^2 (3 r + 8 u)) a e^{-u^2/k^2} / (k^6 r^7)
182 : * \f]
183 : *
184 : * where \f$r\f$ is the extraction radius, \f$u\f$ the retarded time, \f$a\f$
185 : * the amplitude of the pulse and \f$k\f$ the duration of the pulse.
186 : */
187 0 : double dr_pulse_profile_coefficient_b(double time) const;
188 :
189 : /*
190 : * The radial derivative of the C coefficient is calculated as
191 : *
192 : * \f[
193 : * \partial_r C + \partial_t C = -(1/4) (105 * k^8 + 16 r^4 u^4
194 : * - 16 k^2 r^3 u^2 (3 r + 2 u) - 6 k^6 r (9 r + 28 u)
195 : * + 12 k^4 r^2 (r^2 + 4 r u + 9 u^2))
196 : * a e^{-u^2/k^2} / (k^8 r^6)
197 : * \f]
198 : *
199 : * where \f$r\f$ is the extraction radius, \f$u\f$ the retarded time, \f$a\f$
200 : * the amplitude of the pulse and \f$k\f$ the duration of the pulse.
201 : */
202 0 : double dr_pulse_profile_coefficient_c(double time) const;
203 :
204 0 : static DataVector sin_theta(size_t l_max);
205 :
206 0 : static DataVector cos_theta(size_t l_max);
207 :
208 : protected:
209 : /// A no-op as the Teukolsky wave solution does not have substantial
210 : /// shared computation to prepare before the separate component calculations.
211 1 : void prepare_solution(const size_t /*output_l_max*/,
212 : const double /*time*/) const override {}
213 :
214 : /*!
215 : * \brief Compute the spherical coordinate metric from the closed-form
216 : * perturbative Teukolsky wave metric.
217 : *
218 : * \details The specific outgoing wave selected in this analytic solution is
219 : * constructed from a (2,0) mode as in \cite Barkett2019uae, and takes the
220 : * form
221 : *
222 : * \f{align*}{
223 : * g_{tt} &= -1\\
224 : * g_{rr} &= (1 + A f_{rr}) \\
225 : * g_{r \theta} &= 2 B f_{r \theta} r\\
226 : * g_{\theta \theta} &= (1 + C f_{\theta \theta}^{(C)}
227 : * + A f_{\theta \theta}^{(A)}) r^2\\
228 : * g_{\phi \phi} &= (1 + C f_{\phi \phi}^{(C)}
229 : * + A f_{\phi \phi}^{(A)}) r^2 \sin^2 \theta\\
230 : * \f}
231 : *
232 : * and all other components vanish. The angular factors generated by the
233 : * choice of spin-weighted spherical harmonic are
234 : *
235 : * \f{align*}{
236 : * f_{rr} &= 2 - 3 \sin^2 \theta \\
237 : * f_{r \theta} &= -3 \sin \theta \cos \theta \\
238 : * f_{\theta \theta}^{(C)} &= 3 \sin^2 \theta \\
239 : * f_{\theta \theta}^{(A)} &= -1 \\
240 : * f_{\phi \phi}^{(C)} &= - 3 \sin^2 \theta \\
241 : * f_{\phi \phi}^{(A)} &= 3 \sin^2 \theta -1,
242 : * \f}
243 : *
244 : * the radial and time dependent factors are
245 : *
246 : * \f{align*}{
247 : * A &= 3 \left(\frac{\partial_u^2 F(u)}{r^3}
248 : * + \frac{3 \partial_u F(u)}{r^4} + \frac{3 F(u)}{r^5} \right),\\
249 : * B &= - \left(\frac{\partial_u^3 F(u)}{r^2}
250 : * + \frac{3 \partial_u^2 F(u)}{r^3} + \frac{6 \partial_uF(u)}{r^4}
251 : * + \frac{6 F(u)}{r^5}\right), \\
252 : * C &= \frac{1}{4} \left(\frac{\partial_u^4 F(u)}{r}
253 : * + \frac{2 \partial_u^3 F(u)}{r^2} + \frac{9 \partial_u^2 F(u)}{r^3}
254 : * + \frac{21 \partial_u F(u)}{r^4} + \frac{21 F(u)}{r}\right),
255 : * \f}
256 : *
257 : * and the pulse profile is
258 : *
259 : * \f[
260 : * F(u) = a e^{-u^2 /k^2}.
261 : * \f]
262 : *
263 : * So, the pulse profile factors expand to
264 : *
265 : * \f{align*}{
266 : * A &= \frac{3 a e^{-u^2/k^2}}{k^4 r^5} \left(3 k^4 + 4 r^2 u^2
267 : * - 2 k^2 r (r + 3 u)\right),\\
268 : * B &= \frac{2 a e^{-u^2/k^2}}{k^6 r^5} \left(-3 k^6 + 4 r^3 u^3
269 : * - 6 k^2 r^2 u (r + u) + 3 k^4 r (r + 2 u)\right), \\
270 : * C &= \frac{a e^{-u^2/k^2}}{4 k^8 r^5} \left(21 k^8 + 16 r^4 u^4
271 : * - 16 k^2 r^3 u^2 (3 r + u) - 6 k^6 r (3 r + 7 u)
272 : * + 12 k^4 r^2 (r^2 + 2 r u + 3 u^2)\right),
273 : * \f}
274 : *
275 : * \note The \f$\phi\f$ components are returned in a form for which the
276 : * \f$\sin(\theta)\f$ factors are omitted, assuming that derivatives and
277 : * Jacobians will be applied similarly omitting those factors (and therefore
278 : * improving precision of the tensor expression). If you require the
279 : * \f$\sin(\theta)\f$ factors, be sure to put them in by hand in the calling
280 : * code.
281 : */
282 1 : void spherical_metric(
283 : gsl::not_null<
284 : tnsr::aa<DataVector, 3, ::Frame::Spherical<::Frame::Inertial>>*>
285 : spherical_metric,
286 : size_t l_max, double time) const override;
287 :
288 : /*!
289 : * \brief Compute the radial derivative of the spherical coordinate metric
290 : * from the closed-form perturbative Teukolsky wave metric.
291 : *
292 : * \details The specific outgoing wave selected in this analytic solution is
293 : * constructed from a (2,0) mode as in \cite Barkett2019uae, and takes the
294 : * form
295 : *
296 : * \f{align*}{
297 : * \partial_r g_{rr} &= f_{r r} \partial_r A \\
298 : * \partial_r g_{r \theta} &= f_{r \theta} (B + r \partial_r B)\\
299 : * \partial_r g_{\theta \theta} &= 2 (1 + C f_{\theta \theta}^{(C)}
300 : * + A f_{\theta \theta}^{(A)}) r
301 : * + (\partial_r C f_{\theta \theta}^{(C)}
302 : * + \partial_r A f_{\theta \theta}^{(A)}) r^2 \\
303 : * \partial_r g_{\phi \phi} &= 2 (1 + C f_{\phi \phi}^{(C)}
304 : * + A f_{\phi \phi}^{(A)}) r \sin^2 \theta
305 : * + (\partial_r C f_{\phi \phi}^{(C)}
306 : * + \partial_r A f_{\phi \phi}^{(A)}) r^2 \sin^2 \theta\\
307 : * \f}
308 : *
309 : * and all other components vanish. The angular factors \f$f_{a b}\f$ and the
310 : * metric component functions \f$A, B,\f$ and \f$C\f$ are defined as in
311 : * `TeukolskyWave::spherical_metric()`.
312 : * The radial derivatives of the pulse profile functions are obtained by:
313 : *
314 : * \f{align*}{
315 : * \partial_r A + \partial_t A &= \frac{-9 a e^{-u^2/k^2}}{k^4 r^6} \left(
316 : * 5 k^4 + 4 r^2 u^2 - 2 k^2 r (r + 4 u)\right), \\
317 : * \partial_r B + \partial_t B &= \frac{2 a e^{-u^2/k^2}}{k^6 r^6} \left(
318 : * 15 k^6 - 8 r^3 u^3 + 6 k^2 r^2 u (2 r + 3 u)
319 : * - 3 k^4 r (3 r + 8 u)\right), \\
320 : * \partial_r C + \partial_t C &= \frac{-a e^{-u^2/k^2}}{4 k^8 r^6} \left(
321 : * 105 k^8 + 16 k^4 u^4 - 16 k^2 r^3 u^2 (3 r + 2 u) - 6 k^6 r (9 r + 28 u)
322 : * + 12 k^4 r^2 (r^2 + 4 r u + 9 u^2)\right),
323 : * \f}
324 : *
325 : * and the time derivatives of the pulse profile functions are given in
326 : * `TeukolskyWave::dt_spherical_metric()`.
327 : *
328 : * \note The \f$\phi\f$ components are returned in a form for which the
329 : * \f$\sin(\theta)\f$ factors are omitted, assuming that derivatives and
330 : * Jacobians will be applied similarly omitting those factors (and therefore
331 : * improving precision of the tensor expression). If you require the
332 : * \f$\sin(\theta)\f$ factors, be sure to put them in by hand in the calling
333 : * code.
334 : */
335 1 : void dr_spherical_metric(
336 : gsl::not_null<
337 : tnsr::aa<DataVector, 3, ::Frame::Spherical<::Frame::Inertial>>*>
338 : dr_spherical_metric,
339 : size_t l_max, double time) const override;
340 :
341 : /*!
342 : * \brief Compute the time derivative of the spherical coordinate metric
343 : * from the closed-form perturbative Teukolsky wave metric.
344 : *
345 : * \details The specific outgoing wave selected in this analytic solution is
346 : * constructed from a (2,0) mode as in \cite Barkett2019uae, and takes the
347 : * form
348 : *
349 : * \f{align*}{
350 : * \partial_t g_{rr} &= f_{r r} \partial_t A \\
351 : * \partial_t g_{r \theta} &= f_{r \theta} r \partial_t B\\
352 : * \partial_t g_{\theta \theta} &=
353 : * (\partial_t C f_{\theta \theta}^{(C)}
354 : * + \partial_t A f_{\theta \theta}^{(A)}) r^2 \\
355 : * \partial_t g_{\phi \phi} &= (\partial_t C f_{\phi \phi}^{(C)}
356 : * + \partial_t A f_{\phi \phi}^{(A)}) r^2 \sin^2 \theta\\
357 : * \f}
358 : *
359 : * and all other components vanish. The angular factors \f$f_{a b}\f$ and the
360 : * metric component functions \f$A, B,\f$ and \f$C\f$ are defined as in
361 : * `TeukolskyWave::spherical_metric()`.
362 : * The time derivatives of the pulse profile functions are:
363 : *
364 : * \f{align*}{
365 : * \partial_t A &= \frac{-2 u}{k^2} A + \frac{3 a e^{-u^2/k^2}}{k^4 r^5}
366 : * \left( 8 r^2 u - 6 k^2 r \right), \\
367 : * \partial_t B &= \frac{-2 u}{k^2} B + \frac{2 a e^{-u^2/k^2}}{k^6 r^5}
368 : * \left(12 r^3 u^2 - 6 k^2 r^2 (r + 2 u) + 6 k^4 r\right), \\
369 : * \partial_t C &= \frac{-2 u}{k^2} C + \frac{-a e^{-u^2/k^2}}{4 k^8 r^5}
370 : * \left(64 k^4 u^3 - 16 k^2 r^3 u (6 r + 3 u) - 42 k^6 r
371 : * + 12 k^4 r^2 (2 r + 6 u)\right),
372 : * \f}
373 : *
374 : * \note The \f$\phi\f$ components are returned in a form for which the
375 : * \f$\sin(\theta)\f$ factors are omitted, assuming that derivatives and
376 : * Jacobians will be applied similarly omitting those factors (and therefore
377 : * improving precision of the tensor expression). If you require the
378 : * \f$\sin(\theta)\f$ factors, be sure to put them in by hand in the calling
379 : * code.
380 : */
381 1 : void dt_spherical_metric(
382 : gsl::not_null<
383 : tnsr::aa<DataVector, 3, ::Frame::Spherical<::Frame::Inertial>>*>
384 : dt_spherical_metric,
385 : size_t l_max, double time) const override;
386 :
387 0 : using WorldtubeData::variables_impl;
388 :
389 1 : using SphericalMetricData::variables_impl;
390 :
391 : /*!
392 : * \brief Compute the news associated with the (2,0)-mode Teukolsky wave.
393 : *
394 : * \details The value of the news is
395 : *
396 : * \f{align*}{
397 : * N = \frac{3 \sin^2 \theta}{4} \partial_u^5 F(u)
398 : * \f}
399 : *
400 : * where \f$F(u)\f$ is the pulse profile, taken to be
401 : *
402 : * \f[
403 : * F(u) = a e^{-u^2 /k^2},
404 : * \f]
405 : *
406 : * So, the news expands to
407 : *
408 : * \f[
409 : * N = -\frac{6 a e^{-u^2/k^2} u}{k^{10}} \left(15 k^4 - 20 k^2 u^2
410 : * + 4 u^4\right)
411 : * \f]
412 : *
413 : * in this analytic solution.
414 : */
415 1 : void variables_impl(
416 : gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, -2>>*> news,
417 : size_t l_max, double time,
418 : tmpl::type_<Tags::News> /*meta*/) const override;
419 :
420 0 : double amplitude_ = std::numeric_limits<double>::signaling_NaN();
421 0 : double duration_ = std::numeric_limits<double>::signaling_NaN();
422 : };
423 : } // namespace Solutions
424 : } // namespace Cce
|
