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