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