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