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