RobinsonTrautman.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 <memory>
9 #include <vector>
10 
11 #include "DataStructures/ComplexDataVector.hpp"
12 #include "DataStructures/SpinWeighted.hpp"
14 #include "Evolution/Systems/Cce/AnalyticSolutions/SphericalMetricData.hpp"
15 #include "Evolution/Systems/Cce/AnalyticSolutions/WorldtubeData.hpp"
16 #include "NumericalAlgorithms/OdeIntegration/OdeIntegration.hpp"
17 #include "Options/Options.hpp"
19 #include "Time/History.hpp"
21 #include "Utilities/Gsl.hpp"
22 #include "Utilities/Literals.hpp"
23 #include "Utilities/TMPL.hpp"
24 
25 namespace Cce {
26 namespace Solutions {
27 /*!
28  * \brief An analytic solution representing a specialization of the radiative
29  * Robinson-Trautman solution described in \cite Derry1970fk.
30  *
31  * \details This solution is not quite analytic, in the sense that there is a
32  * single scalar field that must be evolved. Ultimately, it is a partial
33  * specialization of the Characteristic equations such that \f$J = 0\f$ and the
34  * evolution equations have been manipulated to give a time evolution equation
35  * for \f$e^{-2 \beta}\f$, which is equivalent to the Robinson-Trautman scalar
36  * \f$\omega_{\text{RT}}\f$ (denoted \f$W\f$ in \cite Derry1970fk -- we deviate
37  * from their notation because the symbol \f$W\f$ is already used elsewhere in
38  * the CCE system).
39  *
40  * \note The value of \f$\omega_{\text{RT}}\f$ should be real and near 1, which
41  * is imposed by the solution itself, any modes specified in the input file are
42  * treated as perturbations to the leading value of 1.0 for the
43  * Robinson-Trautman scalar \f$\omega_{\text{RT}}\f$.
44  */
46  struct InitialModes {
48  static constexpr Options::String help{
49  "The initial modes of the Robinson-Trautman scalar, denoted W in "
50  "[Derry 1970] and omega_RT in the rendered documentation. "
51  "These are taken in ascending l and m order, m varies fastest. Note "
52  "that the modes are treated as perturbations to the leading-order "
53  "solution of 1.0 for omega_RT and only the real part of the field is "
54  "used."};
55  };
57  using type = double;
58  static constexpr Options::String help{
59  "The extraction radius of the spherical solution"};
60  static type lower_bound() noexcept { return 0.0; }
61  static type suggested_value() noexcept { return 20.0; }
62  };
63  struct LMax {
64  using type = size_t;
65  static constexpr Options::String help{
66  "The maximum l value for the internal computation of the analytic "
67  "solution"};
68  static type lower_bound() noexcept { return 4; }
69  };
70  struct Tolerance {
71  using type = double;
72  static constexpr Options::String help{
73  "The tolerance for the time evolution part of the calculation of the "
74  "semi-analytic Robinson-Trautman solution"};
75  static type lower_bound() noexcept { return 0.0; }
76  static type suggested_value() noexcept { return 1.0e-11; }
77  };
78  struct StartTime {
79  using type = double;
80  static constexpr Options::String help{
81  "The starting time for the Robinson-Trautman evolution"};
82  static type lower_bound() noexcept { return 0.0; }
83  static type suggested_value() noexcept { return 0.0; }
84  };
85 
86  using options =
87  tmpl::list<InitialModes, ExtractionRadius, LMax, Tolerance, StartTime>;
88 
89  static constexpr Options::String help = {
90  "Analytic solution representing worldtube data for the nonlinear "
91  "semi-analytic Robinson-Trautman metric, which requires a single "
92  "scalar on the boundary to be evolved to determine the metric"};
93 
94  WRAPPED_PUPable_decl_template(RobinsonTrautman); // NOLINT
95 
96  explicit RobinsonTrautman(CkMigrateMessage* msg) noexcept
97  : SphericalMetricData(msg) {}
98 
99  RobinsonTrautman() = default;
100 
101  RobinsonTrautman(std::vector<std::complex<double>> initial_modes,
102  double extraction_radius, size_t l_max, double tolerance,
103  double start_time, const Options::Context& context);
104 
105  std::unique_ptr<WorldtubeData> get_clone() const noexcept override;
106 
107  void pup(PUP::er& p) noexcept override;
108 
109  bool use_noninertial_news() const noexcept override { return true; }
110 
111  private:
112  void du_rt_scalar(
114  const SpinWeighted<ComplexDataVector, 0>& rt_scalar) const noexcept;
115 
116  void du_bondi_w(
118  const Scalar<SpinWeighted<ComplexDataVector, 0>>& local_du_rt_scalar,
119  const Scalar<SpinWeighted<ComplexDataVector, 0>>& rt_scalar)
120  const noexcept;
121 
122  void bondi_u(
124  const Scalar<SpinWeighted<ComplexDataVector, 0>>& rt_scalar)
125  const noexcept;
126 
127  void bondi_w(
129  const Scalar<SpinWeighted<ComplexDataVector, 0>>& rt_scalar)
130  const noexcept;
131 
132  void dr_bondi_w(
134  const Scalar<SpinWeighted<ComplexDataVector, 0>>& rt_scalar)
135  const noexcept;
136 
137  void initialize_stepper_from_start() noexcept;
138 
139  protected:
140  /*!
141  * \brief The Robinson-Trautman solution performs the time-stepping to advance
142  * the internal member scalar used to generate the metric solution to the
143  * correct state for `time`.
144  *
145  * \details The generating scalar \f$\omega_{\text{RT}}\f$ is evolved
146  * using equation (2.5) from \cite Derry1970fk (manipulated to a form
147  * convenient for our numerical utilities)
148  *
149  * \f[
150  * \partial_u \omega_{\text{RT}} = -
151  * \left(\omega^4_{\text{RT}} \eth^2 \bar \eth^2 \omega_{\text{RT}}
152  * - \omega_{\text{RT}}^3 (\eth^2 \omega_{\text{RT}})
153  * (\bar \eth^2 \omega_{\text{RT}}) \right)
154  * \f]
155  *
156  * As the scalar \f$\omega_{\text{RT}}\f$ is evolved, it is filtered by
157  * zeroing the highest two angular modes.
158  */
159  void prepare_solution(size_t output_l_max,
160  double time) const noexcept override;
161 
162  /*!
163  * \brief Compute the spherical coordinate metric of the Robinson-Trautman
164  * solution generated by the time-evolved scalar \f$\omega_{\text{RT}}\f$.
165  *
166  * \details The spacetime metric of the Robinson-Trautman solution can be
167  * expressed as a specialization of the Bondi-Sachs metric (note the metric
168  * signature change as compared to equation (1.2) from \cite Derry1970fk)
169  *
170  * \f[
171  * ds^2 = -((r W + 1) \omega_{\text{RT}} - r^2 U \bar U) (dt - dr)^2
172  * - 2 \omega_{\text{RT}} (dt - dr) dr
173  * - 2 r^2 U^A q_{AB} dx^B (dt - dr) + r^2 q_{A B} dx^A dx^B,
174  * \f]
175  *
176  * where \f$q_{A B}\f$ represents the angular unit sphere metric, and the
177  * remaining Bondi-Sachs scalars and angular tensors are defined in terms of
178  * the Robinson-Trautman scalar \f$\omega_{\text{RT}}\f$
179  *
180  * \f{align*}{
181  * W &= \frac{1}{r}\left(\omega_{\text{RT}} + \eth \bar \eth
182  * \omega_{\text{RT}} - 1\right) - \frac{2}{r^2 \omega_{\text{RT}}^2}\\
183  * U &\equiv U^A q_A = \frac{\eth \omega_{\text{RT}}}{r}.
184  * \f}
185  * and \f$q_A\f$ is the angular dyad on the unit sphere.
186  *
187  * The angular part of the metric can be expressed in terms of the \f$U\f$
188  * scalar as
189  *
190  * \f{align*}{
191  * g_{u \theta} &= r^2 \Re U\\
192  * g_{u \phi} &= r^2 \Im U
193  * \f}
194  */
195  void spherical_metric(
196  gsl::not_null<
197  tnsr::aa<DataVector, 3, ::Frame::Spherical<::Frame::Inertial>>*>
199  size_t l_max, double time) const noexcept override;
200 
201  /*!
202  * \brief Compute radial derivative of the spherical coordinate metric of the
203  * Robinson-Trautman solution generated by the time-evolved scalar
204  * \f$\omega_{\text{RT}}\f$.
205  *
206  * \details The radial derivative (at constant time t) of the
207  * Robinson-Trautman solution is obtained by differentiating the expressions
208  * from the documentation for `RobinsonTrautman::spherical_metric()`:
209  *
210  * \f{align*}{
211  * (\partial_r g_{a b} + \partial_t g_{a b}) dx^a dx^b =
212  * - (\omega_{\text{RT}} (r \partial_r W + W) - 2 r U \bar U
213  * - r^2 (\bar U\partial_r U + U \partial_r \bar U)) (dt - dr)^2
214  * - 2 r U^A q_{A B} dx^B (dt - dr) + 2 r q_{A B} dx^A dx^B
215  * \f}
216  *
217  * where \f$q_{A B}\f$ represents the angular unit sphere metric, and the
218  * remaining Bondi-Sachs scalars and angular tensors are defined in terms of
219  * the Robinson-Trautman scalar \f$\omega_{\text{RT}}\f$
220  *
221  * \f{align*}{
222  * W &= \frac{1}{r}\left(\omega_{\text{RT}}
223  * + \eth \bar \eth \omega_{\text{RT}} - 1\right)
224  * - \frac{2}{r^2 \omega_{\text{RT}}^2}\\
225  * \partial_r W &= -\frac{1}{r^2} \left(\omega_{\text{RT}}
226  * + \eth \bar \eth \omega_{\text{RT}} - 1\right)
227  * + \frac{4}{r^3 \omega_{\text{RT}}^2}\\
228  * U &\equiv U^A q_A = \frac{\eth \omega_{\text{RT}}}{r}.
229  * \f}
230  * and \f$q_A\f$ is the angular dyad on the unit sphere. The Robinson-Trautman
231  * scalar \f$\omega_{\text{RT}}\f$ is independent of the Bondi radius \f$r\f$,
232  * so all radial derivatives of \f$\omega_{\text{RT}}\f$ have been dropped
233  */
234  void dr_spherical_metric(
235  gsl::not_null<
236  tnsr::aa<DataVector, 3, ::Frame::Spherical<::Frame::Inertial>>*>
238  size_t l_max, double time) const noexcept override;
239 
240  /*!
241  * \brief Compute time derivative of the spherical coordinate metric of the
242  * Robinson-Trautman solution generated by the time-evolved scalar
243  * \f$\omega_{\text{RT}}\f$.
244  *
245  * \details The time derivative of the Robinson-Trautman solution is obtained
246  * by differentiating the expressions from the documentation for
247  * `RobinsonTrautman::spherical_metric()`:
248  *
249  * \f{align*}{
250  * \partial_t g_{a b} dx^a dx^b = -( \partial_u \omega_{\text{RT}} (r W + 1)
251  * + \omega_{\text{RT}} \partial_u W
252  * - r^2 (\bar U \partial_u U + U \partial_u \bar U)) (dt - dr)^2
253  * - 2 \partial_u \omega_{\text{RT}} (dt - dr) dr
254  * - 2 r^2 \partial_u U^A q_{AB} dx^B (dt - dr),
255  * \f}
256  *
257  * where \f$q_{A B}\f$ represents the angular unit sphere metric, and the
258  * remaining Bondi-Sachs scalars and angular tensors are defined in terms of
259  * the Robinson-Trautman scalar \f$\omega_{\text{RT}}\f$
260  *
261  * \f{align*}{
262  * W &= \frac{1}{r}\left(\omega_{\text{RT}}
263  * + \eth \bar \eth \omega_{\text{RT}} - 1\right)
264  * - \frac{2}{r^2 \omega_{\text{RT}}^2}\\
265  * \partial_u W &= \frac{1}{r}\left(\partial_u \omega_{\text{RT}}
266  * + \eth \bar \eth \partial_u \omega_{\text{RT}}\right)
267  * + \frac{4 \partial_u \omega_{\text{RT}}}{r^2 \omega_{\text{RT}}^3} \\
268  * \partial_u U &= q_A \partial_u U^A = \frac{\eth \partial_u
269  * \omega_{\text{RT}}}{r}, \f}
270  *
271  * and \f$q_A\f$ is the angular dyad on the unit sphere; and the time
272  * derivative of the Robinson-Trautman scalar \f$\omega_{\text{RT}}\f$ is
273  *
274  * \f[
275  * \partial_u \omega_{\text{RT}} =
276  * \frac{1}{12} \left(-\omega^4_{\text{RT}} \eth^2 \bar \eth^2
277  * \omega_{\text{RT}} + \omega_{\text{RT}}^3 (\eth^2 \omega_{\text{RT}})
278  * (\bar \eth^2 \omega_{\text{RT}}) \right)
279  * \f]
280  */
281  void dt_spherical_metric(
282  gsl::not_null<
283  tnsr::aa<DataVector, 3, ::Frame::Spherical<::Frame::Inertial>>*>
285  size_t l_max, double time) const noexcept override;
286 
287  /*!
288  * \brief Compute the news associated with the Robinson-Trautman solution
289  * generated by the time-evolved scalar \f$\omega_{\text{RT}}\f$.
290  *
291  * \details The Bondi-Sachs news in the Robinson-Trautman solution is
292  *
293  * \f{align*}{
294  * N = \frac{\bar \eth \bar \eth \omega_{\text{RT}}}{\omega_{\text{RT}}}
295  * \f}
296  */
297  void variables_impl(
298  gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, -2>>*> News,
299  size_t l_max, double time,
300  tmpl::type_<Tags::News> /*meta*/) const noexcept override;
301 
302  using WorldtubeData::variables_impl;
303 
304  using SphericalMetricData::variables_impl;
305  private:
306  mutable std::pair<double, double> step_range_;
307 
308  mutable boost::numeric::odeint::dense_output_runge_kutta<
309  boost::numeric::odeint::controlled_runge_kutta<
310  boost::numeric::odeint::runge_kutta_dopri5<ComplexDataVector>>>
311  stepper_;
312 
313  mutable Scalar<SpinWeighted<ComplexDataVector, 0>> dense_output_rt_scalar_;
314  mutable Scalar<SpinWeighted<ComplexDataVector, 0>> dense_output_du_rt_scalar_;
315  size_t l_max_ = 0;
316  mutable double prepared_time_ = std::numeric_limits<double>::signaling_NaN();
317  double tolerance_ = std::numeric_limits<double>::signaling_NaN();
318  double start_time_ = std::numeric_limits<double>::signaling_NaN();
319  std::vector<std::complex<double>> initial_modes_;
320 };
321 } // namespace Solutions
322 } // namespace Cce
Cce::Solutions::RobinsonTrautman::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 of the Robinson-Trautman solution generated by the time-evolv...
Cce::Solutions::RobinsonTrautman::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 radial derivative of the spherical coordinate metric of the Robinson-Trautman solution genera...
CharmPupable.hpp
Cce::Solutions::RobinsonTrautman::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 Robinson-Trautman solution generated by the time-evolved scalar ...
Cce::Solutions::RobinsonTrautman::StartTime
Definition: RobinsonTrautman.hpp:78
Cce::Solutions::RobinsonTrautman::prepare_solution
void prepare_solution(size_t output_l_max, double time) const noexcept override
The Robinson-Trautman solution performs the time-stepping to advance the internal member scalar used ...
Literals.hpp
Options.hpp
Cce::Solutions::RobinsonTrautman::Tolerance
Definition: RobinsonTrautman.hpp:70
vector
Options::Context
Definition: Options.hpp:41
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::RobinsonTrautman::LMax
Definition: RobinsonTrautman.hpp:63
cstddef
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
memory
tnsr
Type aliases to construct common Tensors.
Definition: TypeAliases.hpp:31
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
Gsl.hpp
Cce::Solutions::RobinsonTrautman::ExtractionRadius
Definition: RobinsonTrautman.hpp:56
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
Cce::Solutions::RobinsonTrautman::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 time derivative of the spherical coordinate metric of the Robinson-Trautman solution generate...
Frame
Definition: IndexType.hpp:36
Tensor.hpp
RungeKutta3.hpp
complex
Cce::Solutions::RobinsonTrautman::InitialModes
Definition: RobinsonTrautman.hpp:46
ComplexDataVector
Stores a collection of complex function values.
Definition: ComplexDataVector.hpp:53
std::unique_ptr
gsl
Implementations from the Guideline Support Library.
Definition: ReadSpecPiecewisePolynomial.hpp:11
Cce::Solutions::RobinsonTrautman
An analytic solution representing a specialization of the radiative Robinson-Trautman solution descri...
Definition: RobinsonTrautman.hpp:45
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13