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