BoundaryData.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 
9 #include "DataStructures/DataBox/PrefixHelpers.hpp"
11 #include "DataStructures/DataVector.hpp"
12 #include "DataStructures/SpinWeighted.hpp"
15 #include "Evolution/Systems/Cce/BoundaryDataTags.hpp"
16 #include "Evolution/Systems/Cce/Tags.hpp"
17 #include "NumericalAlgorithms/Spectral/SwshCollocation.hpp"
18 #include "NumericalAlgorithms/Spectral/SwshDerivatives.hpp"
21 
22 /// \cond
23 class DataVector;
24 class ComplexDataVector;
25 /// \endcond
26 
27 namespace Cce {
28 
29 /*!
30  * \brief Constructs the collocation values for \f$\cos(\phi)\f$,
31  * \f$\cos(\theta)\f$, \f$\sin(\phi)\f$, and \f$\sin(\theta)\f$, returned by
32  * `not_null` pointer in that order.
33  *
34  * \details These are needed for coordinate transformations from the input
35  * Cartesian-like coordinates.
36  */
37 void trigonometric_functions_on_swsh_collocation(
41  gsl::not_null<Scalar<DataVector>*> sin_theta, size_t l_max) noexcept;
42 
43 /*!
44  * \brief Creates both the Jacobian and inverse Jacobian between Cartesian and
45  * spherical coordinates, and the coordinates themselves
46  *
47  * \details The `cartesian_to_spherical_jacobian` is
48  * \f$dx^i/d\tilde{x}^{\tilde j}\f$,
49  * where the Cartesian components are in order \f$x^i = \{x, y, z\}\f$
50  * and the spherical coordinates are
51  * \f$\tilde{x}^{\tilde j} = \{r, \theta, \phi\}\f$.
52  * The Cartesian coordinates given are the standard unit sphere coordinates:
53  *
54  * \f{align*}{
55  * x &= \cos(\phi) \sin(\theta)\\
56  * y &= \sin(\phi) \sin(\theta)\\
57  * z &= \cos(\theta)
58  * \f}
59  *
60  * \note These Jacobians are adjusted to improve regularity near the pole, in
61  * particular the \f$\partial \phi / \partial x^i\f$ components have been scaled
62  * by \f$\sin \theta\f$ (omitting a \f$1/\sin(\theta)\f$) and the
63  * \f$\partial x^i/\partial \phi\f$ components have been scaled by
64  * \f$1/\sin(\theta)\f$ (omitting a \f$\sin(\theta)\f$). The reason is that in
65  * most careful calculations, these problematic sin factors can actually be
66  * omitted because they cancel. In cases where they are actually required, they
67  * must be put in by hand.
68  */
69 void cartesian_to_spherical_coordinates_and_jacobians(
70  gsl::not_null<tnsr::I<DataVector, 3>*> unit_cartesian_coords,
71  gsl::not_null<SphericaliCartesianJ*> cartesian_to_spherical_jacobian,
73  inverse_cartesian_to_spherical_jacobian,
74  const Scalar<DataVector>& cos_phi, const Scalar<DataVector>& cos_theta,
75  const Scalar<DataVector>& sin_phi, const Scalar<DataVector>& sin_theta,
76  double extraction_radius) noexcept;
77 
78 /*
79  * \brief Compute \f$g_{i j}\f$, \f$g^{i j}\f$, \f$\partial_i g_{j k}\f$, and
80  * \f$\partial_t g_{i j}\f$ from input libsharp-compatible modal spatial
81  * metric quantities.
82  *
83  * \details This function interpolates the modes of
84  * input \f$g_{ij}\f$, \f$\partial_r g_{i j}\f$, and \f$\partial_r g_{i j}\f$ to
85  * the libsharp-compatible grid. This function then applies the necessary
86  * jacobian factors and angular derivatives to determine the full \f$\partial_i
87  * g_{j k}\f$.
88  */
89 void cartesian_spatial_metric_and_derivatives_from_modes(
90  gsl::not_null<tnsr::ii<DataVector, 3>*> cartesian_spatial_metric,
91  gsl::not_null<tnsr::II<DataVector, 3>*> inverse_cartesian_spatial_metric,
92  gsl::not_null<tnsr::ijj<DataVector, 3>*> d_cartesian_spatial_metric,
93  gsl::not_null<tnsr::ii<DataVector, 3>*> dt_cartesian_spatial_metric,
95  interpolation_modal_buffer,
97  interpolation_buffer,
99  const tnsr::ii<ComplexModalVector, 3>& spatial_metric_coefficients,
100  const tnsr::ii<ComplexModalVector, 3>& dr_spatial_metric_coefficients,
101  const tnsr::ii<ComplexModalVector, 3>& dt_spatial_metric_coefficients,
102  const CartesianiSphericalJ& inverse_cartesian_to_spherical_jacobian,
103  size_t l_max) noexcept;
104 
105 /*!
106  * \brief Compute \f$\beta^{i}\f$, \f$\partial_i \beta^{j}\f$, and
107  * \f$\partial_t \beta^i\f$ from input libsharp-compatible modal spatial
108  * metric quantities.
109  *
110  * \details This function interpolates the modes of
111  * input \f$\beta^i\f$, \f$\partial_r \beta^i\f$, and \f$\partial_r \beta^i\f$
112  * to the libsharp-compatible grid. This function then applies the necessary
113  * jacobian factors and angular derivatives to determine the full \f$\partial_i
114  * \beta^i\f$.
115  */
116 void cartesian_shift_and_derivatives_from_modes(
117  gsl::not_null<tnsr::I<DataVector, 3>*> cartesian_shift,
118  gsl::not_null<tnsr::iJ<DataVector, 3>*> d_cartesian_shift,
119  gsl::not_null<tnsr::I<DataVector, 3>*> dt_cartesian_shift,
121  interpolation_modal_buffer,
123  interpolation_buffer,
125  const tnsr::I<ComplexModalVector, 3>& shift_coefficients,
126  const tnsr::I<ComplexModalVector, 3>& dr_shift_coefficients,
127  const tnsr::I<ComplexModalVector, 3>& dt_shift_coefficients,
128  const CartesianiSphericalJ& inverse_cartesian_to_spherical_jacobian,
129  size_t l_max) noexcept;
130 
131 /*!
132  * \brief Compute \f$\alpha\f$, \f$\partial_i \alpha\f$, and
133  * \f$\partial_t \beta^i\f$ from input libsharp-compatible modal spatial
134  * metric quantities.
135  *
136  * \details This function interpolates the modes of input \f$\alpha\f$,
137  * \f$\partial_r \alpha\f$, and \f$\partial_r \alpha\f$ to the
138  * libsharp-compatible grid. This function then applies the necessary jacobian
139  * factors and angular derivatives to determine the full \f$\partial_i
140  * \alpha\f$.
141  */
142 void cartesian_lapse_and_derivatives_from_modes(
143  gsl::not_null<Scalar<DataVector>*> cartesian_lapse,
144  gsl::not_null<tnsr::i<DataVector, 3>*> d_cartesian_lapse,
145  gsl::not_null<Scalar<DataVector>*> dt_cartesian_lapse,
147  interpolation_modal_buffer,
149  interpolation_buffer,
151  const Scalar<ComplexModalVector>& lapse_coefficients,
152  const Scalar<ComplexModalVector>& dr_lapse_coefficients,
153  const Scalar<ComplexModalVector>& dt_lapse_coefficients,
154  const CartesianiSphericalJ& inverse_cartesian_to_spherical_jacobian,
155  size_t l_max) noexcept;
156 
157 /*!
158  * \brief Computes the spacetime metric and its first derivative in the
159  * intermediate radial null coordinates
160  *
161  * \details These components are obtained by the steps in
162  * Section II-A of \cite Barkett2019uae, which is based on the computation from
163  * Section 4.3 of \cite Bishop1998uk. The most direct comparison is to be made
164  * with equation (31) of \cite Barkett2019uae, which gives the null metric
165  * components explicitly. The time derivative is then (using notation from
166  * equation (31) of \cite Barkett2019uae):
167  *
168  * \f{align}{
169  * \partial_{\bar u} g_{\bar u \bar \lambda} =
170  * \partial_{\bar u} g_{\bar \lambda \bar \lambda} =
171  * \partial_{\bar u} g_{\bar \lambda \bar A} &= 0 \\
172  * \partial_{\bar u} g_{\bar u \bar u} &=
173  * \partial_{\breve t} g_{\breve t \breve t} \\
174  * \partial_{\bar u} g_{\bar u \bar A} &=
175  * \frac{\partial \breve x^{\breve i}}{\partial \bar x^{\bar A}}\\
176  * g_{\breve i \breve t}
177  * \partial_{\bar u} g_{\bar A \bar B}
178  * &= \frac{\partial \breve x^{\breve i}}{\partial \bar x^{\bar A}}
179  * \frac{\partial \breve x^{\breve j}}{\partial \bar x^{\bar B}}
180  * g_{\breve i \breve j}
181  * \f}
182  */
183 void null_metric_and_derivative(
184  gsl::not_null<tnsr::aa<DataVector, 3, Frame::RadialNull>*> du_null_metric,
185  gsl::not_null<tnsr::aa<DataVector, 3, Frame::RadialNull>*> null_metric,
186  const SphericaliCartesianJ& cartesian_to_spherical_jacobian,
187  const tnsr::aa<DataVector, 3>& dt_spacetime_metric,
188  const tnsr::aa<DataVector, 3>& spacetime_metric) noexcept;
189 
190 /*!
191  * \brief Computes the spatial unit normal vector \f$s^i\f$ to the spherical
192  * worldtube surface and its first time derivative.
193  *
194  * \details Refer to equation (20) of \cite Barkett2019uae for the expression of
195  * the spatial unit normal vector, and equation (23) of \cite Barkett2019uae for
196  * the first time derivative. Refer to \cite Bishop1998uk for more exposition
197  * about the overall construction of the coordinate transformations used for the
198  * intermediate null coordinates.
199  */
200 void worldtube_normal_and_derivatives(
201  gsl::not_null<tnsr::I<DataVector, 3>*> worldtube_normal,
202  gsl::not_null<tnsr::I<DataVector, 3>*> dt_worldtube_normal,
203  const Scalar<DataVector>& cos_phi, const Scalar<DataVector>& cos_theta,
204  const tnsr::aa<DataVector, 3>& spacetime_metric,
205  const tnsr::aa<DataVector, 3>& dt_spacetime_metric,
206  const Scalar<DataVector>& sin_phi, const Scalar<DataVector>& sin_theta,
207  const tnsr::II<DataVector, 3>& inverse_spatial_metric) noexcept;
208 
209 /*!
210  * \brief Computes the null 4-vector \f$l^\mu\f$ on the worldtube surface that
211  * is to be used as the CCE hypersurface generator, and the first time
212  * derivative \f$\partial_u l^\mu\f$.
213  *
214  * \details For mathematical description of our choice of the null generator,
215  * refer to equation (22) of \cite Barkett2019uae, and for the first time
216  * derivative see (25) of \cite Barkett2019uae. Refer to \cite Bishop1998uk for
217  * more exposition about the overall construction of the coordinate
218  * transformations used for the intermediate null coordinates.
219  */
220 void null_vector_l_and_derivatives(
221  gsl::not_null<tnsr::A<DataVector, 3>*> du_null_l,
222  gsl::not_null<tnsr::A<DataVector, 3>*> null_l,
223  const tnsr::I<DataVector, 3>& dt_worldtube_normal,
224  const Scalar<DataVector>& dt_lapse,
225  const tnsr::aa<DataVector, 3>& dt_spacetime_metric,
226  const tnsr::I<DataVector, 3>& dt_shift, const Scalar<DataVector>& lapse,
227  const tnsr::aa<DataVector, 3>& spacetime_metric,
228  const tnsr::I<DataVector, 3>& shift,
229  const tnsr::I<DataVector, 3>& worldtube_normal) noexcept;
230 
231 /*!
232  * \brief Computes the partial derivative of the spacetime metric and inverse
233  * spacetime metric in the intermediate null radial coordinates with respect to
234  * the null generator \f$l^\mu\f$
235  *
236  * \details For full expressions of the \f$l^\mu \partial_\mu g_{a b}\f$ and
237  * \f$l^\mu \partial_\mu g^{a b}\f$ computed in this function, see equation (31)
238  * and (32) of \cite Barkett2019uae. Refer to \cite Bishop1998uk for more
239  * exposition about the overall construction of the coordinate transformations
240  * used for the intermediate null coordinates.
241  */
242 void dlambda_null_metric_and_inverse(
243  gsl::not_null<tnsr::aa<DataVector, 3, Frame::RadialNull>*>
244  dlambda_null_metric,
245  gsl::not_null<tnsr::AA<DataVector, 3, Frame::RadialNull>*>
246  dlambda_inverse_null_metric,
247  const AngulariCartesianA& angular_d_null_l,
248  const SphericaliCartesianJ& cartesian_to_spherical_jacobian,
249  const tnsr::iaa<DataVector, 3>& phi,
250  const tnsr::aa<DataVector, 3>& dt_spacetime_metric,
251  const tnsr::A<DataVector, 3>& du_null_l,
252  const tnsr::AA<DataVector, 3, Frame::RadialNull>& inverse_null_metric,
253  const tnsr::A<DataVector, 3>& null_l,
254  const tnsr::aa<DataVector, 3>& spacetime_metric) noexcept;
255 
256 /*!
257  * \brief Computes the Bondi radius of the worldtube.
258  *
259  * \details Note that unlike the Cauchy coordinate radius, the Bondi radius is
260  * not constant over the worldtube. Instead, it is obtained by the determinant
261  * of the angular part of the metric in the intermediate null coordinates (see
262  * \cite Barkett2019uae).
263  *
264  * \f[
265  * r = \left(\frac{\det g_{A B}}{ q_{A B}}\right)^{1/4},
266  * \f]
267  *
268  * where \f$q_{A B}\f$ is the unit sphere metric.
269  */
270 void bondi_r(
272  const tnsr::aa<DataVector, 3, Frame::RadialNull>& null_metric) noexcept;
273 
274 /*!
275  * \brief Computes the full 4-dimensional partial of the Bondi radius with
276  * respect to the intermediate null coordinates.
277  *
278  * \details The expression evaluated is obtained from differentiating the
279  * determinant equation for `bondi_r`, from (35) of \cite Barkett2019uae :
280  *
281  * \f[
282  * \partial_\alpha r = \frac{r}{4} \left(g^{A B} \partial_\alpha g_{A B}
283  * - \frac{\partial_\alpha \det q_{A B}}{\det q_{A B}}\right)
284  * \f]
285  *
286  * Note that for the angular derivatives, we just numerically differentiate
287  * using the utilities in `Spectral::Swsh::angular_derivative()`. For the time
288  * and radial derivatives, the second term in the above equation vanishes.
289  */
290 void d_bondi_r(
291  gsl::not_null<tnsr::a<DataVector, 3, Frame::RadialNull>*> d_bondi_r,
293  const tnsr::aa<DataVector, 3, Frame::RadialNull>& dlambda_null_metric,
294  const tnsr::aa<DataVector, 3, Frame::RadialNull>& du_null_metric,
295  const tnsr::AA<DataVector, 3, Frame::RadialNull>& inverse_null_metric,
296  size_t l_max) noexcept;
297 
298 /*!
299  * \brief Compute the complex angular dyads used to define the spin-weighted
300  * scalars in the CCE system.
301  *
302  * \details We use the typically chosen angular dyads in CCE
303  * \cite Barkett2019uae \cite Bishop1997ik :
304  *
305  * \f{align*}{
306  * q_A &= \{-1, -i \sin(\theta)\}\\
307  * q^A &= \left\{-1, -i \frac{1}{\sin \theta}\right\}
308  * \f}
309  *
310  * However, to maintain regularity and for compatibility with the more regular
311  * Jacobians from `Cce::cartesian_to_spherical_coordinates_and_jacobians()`, in
312  * the code we omit the factors of \f$\sin \theta\f$ from the above equations.
313  */
314 void dyads(
315  gsl::not_null<tnsr::i<ComplexDataVector, 2, Frame::RadialNull>*> down_dyad,
316  gsl::not_null<tnsr::I<ComplexDataVector, 2, Frame::RadialNull>*>
317  up_dyad) noexcept;
318 
319 /*!
320  * \brief Compute the \f$\beta\f$ (lapse) function for the CCE Bondi-like
321  * metric.
322  *
323  * \details The Bondi-like metric has \f$g^{u r} = - e^{2 \beta}\f$, and the
324  * value of \f$\beta\f$ is obtained from the intermediate null metric by (see
325  * equation (51) of \cite Barkett2019uae) using:
326  *
327  * \f[
328  * \beta = -\frac{1}{2} \ln \partial_{\lambda} r
329  * \f]
330  */
331 void beta_worldtube_data(
333  const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r) noexcept;
334 
335 /*!
336  * \brief Compute the \f$U\f$ (shift) function for the CCE Bondi-like metric.
337  *
338  * \details The Bondi-like metric has \f$g^{r A} = -e^{-2 \beta} U^A\f$, and the
339  * spin-weighted vector \f$U = U^A q_A\f$. The value of \f$U^A\f$ can be
340  * computed from the intermediate null metric quantities (see equation (54) of
341  * \cite Barkett2019uae) using:
342  *
343  * \f[
344  * U = -(\partial_\lambda r g^{\lambda A} + \partial_B r g^{A B}) q_A
345  * / \partial_\lambda r \f]
346  *
347  */
348 void bondi_u_worldtube_data(
350  const tnsr::i<ComplexDataVector, 2, Frame::RadialNull>& dyad,
351  const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
352  const tnsr::AA<DataVector, 3, Frame::RadialNull>&
353  inverse_null_metric) noexcept;
354 
355 /*!
356  * \brief Compute the \f$W\f$ (mass aspect) function for the CCE Bondi-like
357  * metric.
358  *
359  * \details The Bondi-like metric has \f$g^{rr} = e^{-2 \beta}(1 + r W)\f$. The
360  * value of \f$W\f$ can be computed from the null metric quantities (see
361  * equation (55) of \cite Barkett2019uae) using:
362  *
363  * \f[
364  * W = \frac{1}{r} \left(-1
365  * + \frac{g^{\lambda \lambda} (\partial_\lambda r)^2
366  * + 2 \partial_\lambda r \left(\partial_A r g^{\lambda A}
367  * - \partial_u r\right) + \partial_A r \partial_B r g^{A B}}
368  * {\partial_\lambda r}\right) \f]
369  */
370 void bondi_w_worldtube_data(
372  const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
373  const tnsr::AA<DataVector, 3, Frame::RadialNull>& inverse_null_metric,
374  const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r) noexcept;
375 
376 /*!
377  * \brief Compute the \f$J\f$ (intuitively similar to the transverse-traceless
378  * part of the angular metric) function for the CCE Bondi-like metric.
379  *
380  * \details The Bondi-like metric has \f$J = \frac{1}{2 r^2} q^A q^B g_{A B}\f$.
381  * This expression holds both for the right-hand side in the Bondi coordinates
382  * and for the right-hand side in the intermediate null coordinates (see
383  * equation (45) of \cite Barkett2019uae).
384  */
385 void bondi_j_worldtube_data(
387  const tnsr::aa<DataVector, 3, Frame::RadialNull>& null_metric,
389  const tnsr::I<ComplexDataVector, 2, Frame::RadialNull>& dyad) noexcept;
390 
391 /*!
392  * \brief Compute the radial derivative of the angular metric spin-weighted
393  * scalar \f$\partial_r J\f$ in the CCE Bondi-like metric.
394  *
395  * \details The radial derivative of the angular spin-weighted scalar \f$J\f$
396  * can be computed from the null metric components by (c.f. equation (47) of
397  * \cite Barkett2019uae):
398  *
399  * \f[
400  * \partial_r J = \frac{\partial_\lambda J}{\partial_\lambda r} =
401  * \frac{q^A q^B \partial_\lambda g_{A B} / (2 r^2)
402  * - 2 \partial_\lambda r J / r}{\partial_\lambda r}
403  * \f]
404  */
405 void dr_bondi_j(
408  denominator_buffer,
409  const tnsr::aa<DataVector, 3, Frame::RadialNull>& dlambda_null_metric,
410  const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
413  const tnsr::I<ComplexDataVector, 2, Frame::RadialNull>& dyad) noexcept;
414 
415 /*!
416  * \brief Compute the second derivative of the Bondi radius with respect to the
417  * intermediate null coordinate radius \f$\partial_\lambda^2 r\f$.
418  *
419  * \details To determine this second derivative quantity without resorting to
420  * depending on second-derivative metric inputs, we need to take advantage of
421  * one of the Einstein field equations. Combining equations (53) and (52) of
422  * \cite Barkett2019uae, we have:
423  *
424  * \f[
425  * \partial_\lambda^2 r = \frac{-r}{4} \left(
426  * \partial_\lambda J \partial_\lambda \bar J - (\partial_\lambda K)^2\right)
427  * \f],
428  *
429  * where the first derivative of \f$K\f$ can be obtained from \f$K = \sqrt{1 + J
430  * \bar J}\f$ and the first derivative of \f$J\f$ can be obtained from (47) of
431  * \cite Barkett2019uae
432  */
433 void d2lambda_bondi_r(
434  gsl::not_null<Scalar<DataVector>*> d2lambda_bondi_r,
435  const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
436  const Scalar<SpinWeighted<ComplexDataVector, 2>>& dr_bondi_j,
438  const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r) noexcept;
439 
440 /*!
441  * \brief Compute the Bondi metric contribution \f$Q\f$ (radial derivative of
442  * shift).
443  *
444  * \details The definition of \f$Q\f$ in terms of the Bondi metric components is
445  *
446  * \f[
447  * Q = q^A e^{-2 \beta} g_{A B} \partial_r U^B.
448  * \f]
449  *
450  * $Q$ can be derived from the intermediate null metric quantities via (see
451  * equations (56) and (57) of \cite Barkett2019uae)
452  *
453  * \f[
454  * \partial_\lambda U = - \left(\partial_\lambda g^{\lambda A}
455  * + \frac{\partial_A \partial_\lambda r}{\partial_\lambda r} g^{A B}
456  * + \frac{\partial_B r}{\partial_\lambda r} \partial_\lambda g^{A B}\right) q_A
457  * + 2 \partial_\lambda \beta (U + g^{\lambda A} q_A)
458  * \f]
459  *
460  * and
461  *
462  * \f[
463  * Q = r^2 (J \partial_\lambda \bar U + K \partial_\lambda U)
464  * \f]
465  *
466  * also provided is \f$\partial_r U\f$, which is separately useful to cache for
467  * other intermediate steps in the CCE computation.
468  */
469 void bondi_q_worldtube_data(
472  const Scalar<DataVector>& d2lambda_r,
473  const tnsr::AA<DataVector, 3, Frame::RadialNull>&
474  dlambda_inverse_null_metric,
475  const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
476  const tnsr::i<ComplexDataVector, 2, Frame::RadialNull>& dyad,
477  const tnsr::i<DataVector, 2, Frame::RadialNull>& angular_d_dlambda_r,
478  const tnsr::AA<DataVector, 3, Frame::RadialNull>& inverse_null_metric,
481  const Scalar<SpinWeighted<ComplexDataVector, 1>>& bondi_u) noexcept;
482 
483 /*!
484  * \brief Compute the Bondi metric contribution \f$(\partial_u J)_{y} \equiv
485  * H\f$ (the retarded time derivative evaluated at fixed $y$ coordinate) on the
486  * worldtube boundary.
487  *
488  * \details The numerical time derivative (along the worldtube, rather than
489  * along the surface of constant Bondi \f$r\f$) is computed by (see equation
490  * (48) of \cite Barkett2019uae)
491  *
492  * \f[
493  * (\partial_u J)_y = \frac{1}{2 r^2} q^A q^B \partial_u g_{A B}
494  * - \frac{2 \partial_u r}{r} J
495  * \f]
496  *
497  * \note There is the regrettable notation difference with the primary reference
498  * for these formulas \cite Barkett2019uae in that we denote with \f$H\f$ the
499  * time derivative at constant numerical radius, where \cite Barkett2019uae uses
500  * \f$H\f$ to denote the time derivative at constant Bondi radius.
501  */
502 void bondi_h_worldtube_data(
504  const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
506  const tnsr::aa<DataVector, 3, Frame::RadialNull>& du_null_metric,
508  const tnsr::I<ComplexDataVector, 2, Frame::RadialNull>& dyad) noexcept;
509 
510 /*!
511  * \brief Compute the Bondi metric contribution \f$(\partial_u J)_r\f$ (the
512  * retarded time derivative at fixed coordinate $r$) on the worldtube boundary.
513  *
514  * \details The numerical time derivative (along the surface of constant r, not
515  * along the worldtube) is computed by (see equation (50) of
516  * \cite Barkett2019uae)
517  *
518  * \f[
519  * \partial_u J = \frac{1}{2 r^2} q^A q^B \left(\partial_u g_{A B}
520  * - \frac{ \partial_u r}{ \partial_\lambda r} \partial_\lambda g_{A B}\right)
521  * \f]
522  *
523  * \note There is the regrettable notation difference with the primary reference
524  * for these formulas \cite Barkett2019uae in that we denote with \f$H\f$ the
525  * time derivative at constant numerical radius, where \cite Barkett2019uae uses
526  * \f$H\f$ to denote the time derivative at constant Bondi radius.
527  */
528 void du_j_worldtube_data(
530  const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
532  const tnsr::aa<DataVector, 3, Frame::RadialNull>& du_null_metric,
533  const tnsr::aa<DataVector, 3, Frame::RadialNull>& dlambda_null_metric,
535  const tnsr::I<ComplexDataVector, 2, Frame::RadialNull>& dyad) noexcept;
536 
537 namespace Tags {
538 /// The collection of tags mutated by `create_bondi_boundary_data`
539 template <template <typename> class BoundaryPrefix>
541  BoundaryPrefix,
542  tmpl::list<Tags::BondiBeta, Tags::BondiU, Tags::Dr<Tags::BondiU>,
546 } // namespace Tags
547 
548 namespace detail {
549 // the common step between the modal input and the Generalized harmonic input
550 // that performs the final gauge processing to Bondi scalars and places them in
551 // the DataBox.
552 template <typename DataBoxTagList, typename BufferTagList,
553  typename ComplexBufferTagList>
555  const gsl::not_null<db::DataBox<DataBoxTagList>*> bondi_boundary_data,
556  const gsl::not_null<Variables<BufferTagList>*> computation_variables,
557  const gsl::not_null<Variables<ComplexBufferTagList>*> derivative_buffers,
558  const tnsr::aa<DataVector, 3>& dt_spacetime_metric,
559  const tnsr::iaa<DataVector, 3>& phi,
560  const tnsr::aa<DataVector, 3>& spacetime_metric,
561  const tnsr::A<DataVector, 3>& null_l,
562  const tnsr::A<DataVector, 3>& du_null_l,
563  const SphericaliCartesianJ& cartesian_to_spherical_jacobian,
564  const size_t l_max, const double extraction_radius) noexcept {
565  const size_t size = Spectral::Swsh::number_of_swsh_collocation_points(l_max);
566 
567  // unfortunately, because the dyads are not themselves spin-weighted, they
568  // need a separate Variables
569  Variables<tmpl::list<Tags::detail::DownDyad, Tags::detail::UpDyad>>
570  dyad_variables{size};
571 
572  auto& null_metric =
573  get<gr::Tags::SpacetimeMetric<3, Frame::RadialNull, DataVector>>(
574  *computation_variables);
575  auto& du_null_metric = get<
577  *computation_variables);
578  null_metric_and_derivative(
579  make_not_null(&du_null_metric), make_not_null(&null_metric),
580  cartesian_to_spherical_jacobian, dt_spacetime_metric, spacetime_metric);
581 
582  auto& inverse_null_metric =
583  get<gr::Tags::InverseSpacetimeMetric<3, Frame::RadialNull, DataVector>>(
584  *computation_variables);
585 
586  // the below scaling process is used to reduce accumulation of numerical
587  // error in the determinant evaluation
588 
589  // buffer reuse because the scaled null metric is only needed until the
590  // `determinant_and_inverse` call
591  auto& scaled_null_metric =
592  get<gr::Tags::InverseSpacetimeMetric<3, Frame::RadialNull, DataVector>>(
593  *computation_variables);
594  for (size_t i = 0; i < 4; ++i) {
595  for (size_t j = i; j < 4; ++j) {
596  if (i > 1 and j > 1) {
597  scaled_null_metric.get(i, j) =
598  null_metric.get(i, j) / square(extraction_radius);
599  } else if (i > 1 or j > 1) {
600  scaled_null_metric.get(i, j) =
601  null_metric.get(i, j) / extraction_radius;
602  } else {
603  scaled_null_metric.get(i, j) = null_metric.get(i, j);
604  }
605  }
606  }
607  // Allocation
608  const auto scaled_inverse_null_metric =
609  determinant_and_inverse(scaled_null_metric).second;
610  for (size_t i = 0; i < 4; ++i) {
611  for (size_t j = i; j < 4; ++j) {
612  if (i > 1 and j > 1) {
613  inverse_null_metric.get(i, j) =
614  scaled_inverse_null_metric.get(i, j) / square(extraction_radius);
615  } else if (i > 1 or j > 1) {
616  inverse_null_metric.get(i, j) =
617  scaled_inverse_null_metric.get(i, j) / extraction_radius;
618  } else {
619  inverse_null_metric.get(i, j) = scaled_inverse_null_metric.get(i, j);
620  }
621  }
622  }
623 
624  auto& angular_d_null_l =
625  get<Tags::detail::AngularDNullL>(*computation_variables);
626  auto& buffer_for_derivatives =
627  get(get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
629  *derivative_buffers));
630  auto& eth_buffer =
631  get(get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
633  *derivative_buffers));
634  for (size_t a = 0; a < 4; ++a) {
635  buffer_for_derivatives.data() =
636  std::complex<double>(1.0, 0.0) * null_l.get(a);
637  Spectral::Swsh::angular_derivatives<tmpl::list<Spectral::Swsh::Tags::Eth>>(
638  l_max, 1, make_not_null(&eth_buffer), buffer_for_derivatives);
639  angular_d_null_l.get(0, a) = -real(eth_buffer.data());
640  angular_d_null_l.get(1, a) = -imag(eth_buffer.data());
641  }
642 
643  auto& dlambda_null_metric = get<Tags::detail::DLambda<
645  *computation_variables);
646  auto& dlambda_inverse_null_metric = get<Tags::detail::DLambda<
648  *computation_variables);
649  dlambda_null_metric_and_inverse(
650  make_not_null(&dlambda_null_metric),
651  make_not_null(&dlambda_inverse_null_metric), angular_d_null_l,
652  cartesian_to_spherical_jacobian, phi, dt_spacetime_metric, du_null_l,
653  inverse_null_metric, null_l, spacetime_metric);
654 
655  auto& r = db::get<Tags::BoundaryValue<Tags::BondiR>>(*bondi_boundary_data);
656  db::mutate<Tags::BoundaryValue<Tags::BondiR>>(bondi_boundary_data, bondi_r,
657  null_metric);
658 
659  auto& d_r =
660  get<::Tags::spacetime_deriv<Tags::detail::RealBondiR, tmpl::size_t<3>,
661  Frame::RadialNull>>(*computation_variables);
662  d_bondi_r(make_not_null(&d_r), r, dlambda_null_metric, du_null_metric,
663  inverse_null_metric, l_max);
664  db::mutate<Tags::BoundaryValue<Tags::DuRDividedByR>,
666  bondi_boundary_data,
667  [&d_r, &
669  du_r_divided_by_r,
671  du_r) noexcept {
672  get(*du_r_divided_by_r).data() =
673  std::complex<double>{1.0, 0.0} * get<0>(d_r) / get(r).data();
674  get(*du_r).data() = std::complex<double>{1.0, 0.0} * get<0>(d_r);
675  });
676 
677  auto& down_dyad = get<Tags::detail::DownDyad>(dyad_variables);
678  auto& up_dyad = get<Tags::detail::UpDyad>(dyad_variables);
679  dyads(make_not_null(&down_dyad), make_not_null(&up_dyad));
680 
681  db::mutate<Tags::BoundaryValue<Tags::BondiBeta>>(bondi_boundary_data,
682  beta_worldtube_data, d_r);
683 
684  auto& bondi_u =
685  db::get<Tags::BoundaryValue<Tags::BondiU>>(*bondi_boundary_data);
686  db::mutate<Tags::BoundaryValue<Tags::BondiU>>(
687  bondi_boundary_data, bondi_u_worldtube_data, down_dyad, d_r,
688  inverse_null_metric);
689 
690  db::mutate<Tags::BoundaryValue<Tags::BondiW>>(
691  bondi_boundary_data, bondi_w_worldtube_data, d_r, inverse_null_metric, r);
692 
693  auto& bondi_j =
694  db::get<Tags::BoundaryValue<Tags::BondiJ>>(*bondi_boundary_data);
695  db::mutate<Tags::BoundaryValue<Tags::BondiJ>>(
696  bondi_boundary_data, bondi_j_worldtube_data, null_metric, r, up_dyad);
697 
698  auto& dr_j = db::get<Tags::BoundaryValue<Tags::Dr<Tags::BondiJ>>>(
699  *bondi_boundary_data);
700  auto& denominator_buffer =
701  get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
703  *derivative_buffers);
704  db::mutate<Tags::BoundaryValue<Tags::Dr<Tags::BondiJ>>>(
705  bondi_boundary_data, dr_bondi_j, make_not_null(&denominator_buffer),
706  dlambda_null_metric, d_r, bondi_j, r, up_dyad);
707 
708  auto& d2lambda_r = get<
709  Tags::detail::DLambda<Tags::detail::DLambda<Tags::detail::RealBondiR>>>(
710  *computation_variables);
711  d2lambda_bondi_r(make_not_null(&d2lambda_r), d_r, dr_j, bondi_j, r);
712 
713  auto& angular_d_dlambda_r =
714  get<::Tags::deriv<Tags::detail::DLambda<Tags::detail::RealBondiR>,
715  tmpl::size_t<2>, Frame::RadialNull>>(
716  *computation_variables);
717  buffer_for_derivatives.data() = std::complex<double>(1.0, 0.0) * get<1>(d_r);
718  Spectral::Swsh::angular_derivatives<tmpl::list<Spectral::Swsh::Tags::Eth>>(
719  l_max, 1, make_not_null(&eth_buffer), buffer_for_derivatives);
720  angular_d_dlambda_r.get(0) = -real(eth_buffer.data());
721  angular_d_dlambda_r.get(1) = -imag(eth_buffer.data());
722 
723  db::mutate<Tags::BoundaryValue<Tags::BondiQ>,
725  bondi_boundary_data, bondi_q_worldtube_data, d2lambda_r,
726  dlambda_inverse_null_metric, d_r, down_dyad, angular_d_dlambda_r,
727  inverse_null_metric, bondi_j, r, bondi_u);
728 
729  db::mutate<Tags::BoundaryValue<Tags::BondiH>>(
730  bondi_boundary_data, bondi_h_worldtube_data, d_r, bondi_j, du_null_metric,
731  r, up_dyad);
732 
733  db::mutate<Tags::BoundaryValue<Tags::Du<Tags::BondiJ>>>(
734  bondi_boundary_data, du_j_worldtube_data, d_r, bondi_j, du_null_metric,
735  dlambda_null_metric, r, up_dyad);
736 }
737 } // namespace detail
738 
739 /*!
740  * \brief Process the worldtube data from generalized harmonic quantities
741  * to desired Bondi quantities, placing the result in the passed
742  * \ref DataBoxGroup.
743  *
744  * \details
745  * The mathematics are a bit complicated for all of the coordinate
746  * transformations that are necessary to obtain the Bondi gauge quantities.
747  * For full mathematical details, see the documentation for functions in
748  * `BoundaryData.hpp` and \cite Barkett2019uae \cite Bishop1998uk.
749  *
750  * This function takes as input the full set of Generalized harmonic metric data
751  * on a two-dimensional surface of constant \f$r\f$ and \f$t\f$ in numerical
752  * coordinates.
753  *
754  * Sufficient tags to provide full worldtube boundary data at a particular
755  * time are set in `bondi_boundary_data`. In particular, the set of tags in
756  * `Tags::characteristic_worldtube_boundary_tags` in the provided \ref
757  * DataBoxGroup are assigned to the worldtube boundary values associated with
758  * the input metric components.
759  *
760  * The majority of the mathematical transformations are implemented as a set of
761  * individual cascaded functions below. The details of the manipulations that
762  * are performed to the input data may be found in the individual functions
763  * themselves, which are called in the following order:
764  * - `trigonometric_functions_on_swsh_collocation()`
765  * - `gr::shift()`
766  * - `gr::lapse()`
767  * - `worldtube_normal_and_derivatives()`
768  * - `gr::spacetime_normal_vector()`
769  * - `GeneralizedHarmonic::time_deriv_of_lapse()`
770  * - `GeneralizedHarmonic::time_deriv_of_shift()`
771  * - `null_vector_l_and_derivatives()`
772  * - `cartesian_to_spherical_coordinates_and_jacobians()`
773  * - `null_metric_and_derivative()`
774  * - `dlambda_null_metric_and_inverse()`
775  * - `bondi_r()`
776  * - `d_bondi_r()`
777  * - `dyads()`
778  * - `beta_worldtube_data()`
779  * - `bondi_u_worldtube_data()`
780  * - `bondi_w_worldtube_data()`
781  * - `bondi_j_worldtube_data()`
782  * - `dr_bondi_j()`
783  * - `d2lambda_bondi_r()`
784  * - `bondi_q_worldtube_data()`
785  * - `bondi_h_worldtube_data()`
786  * - `du_j_worldtube_data()`
787  */
788 template <typename DataBoxTagList>
790  const gsl::not_null<db::DataBox<DataBoxTagList>*> bondi_boundary_data,
791  const tnsr::iaa<DataVector, 3>& phi, const tnsr::aa<DataVector, 3>& pi,
792  const tnsr::aa<DataVector, 3>& spacetime_metric,
793  const double extraction_radius, const size_t l_max) noexcept {
794  const size_t size = Spectral::Swsh::number_of_swsh_collocation_points(l_max);
795  // Most allocations required for the full boundary computation are merged into
796  // a single, large Variables allocation. There remain a handful of cases in
797  // the computational functions called where an intermediate quantity that is
798  // not re-used is allocated rather than taking a buffer. These cases are
799  // marked with code comments 'Allocation'; In the future, if allocations are
800  // identified as a point to optimize, those buffers may be allocated here and
801  // passed as function arguments
802  Variables<tmpl::list<
803  Tags::detail::CosPhi, Tags::detail::CosTheta, Tags::detail::SinPhi,
804  Tags::detail::SinTheta, Tags::detail::CartesianCoordinates,
805  Tags::detail::CartesianToSphericalJacobian,
806  Tags::detail::InverseCartesianToSphericalJacobian,
813  Tags::detail::WorldtubeNormal, ::Tags::dt<Tags::detail::WorldtubeNormal>,
815  Tags::detail::NullL, ::Tags::dt<Tags::detail::NullL>,
816  // for the detail function called at the end
820  Tags::detail::AngularDNullL,
821  Tags::detail::DLambda<
822  gr::Tags::SpacetimeMetric<3, Frame::RadialNull, DataVector>>,
823  Tags::detail::DLambda<
824  gr::Tags::InverseSpacetimeMetric<3, Frame::RadialNull, DataVector>>,
827  Tags::detail::DLambda<Tags::detail::DLambda<Tags::detail::RealBondiR>>,
829  tmpl::size_t<2>, Frame::RadialNull>>>
830  computation_variables{size};
831 
832  Variables<
833  tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
837  derivative_buffers{size};
838 
839  auto& cos_phi = get<Tags::detail::CosPhi>(computation_variables);
840  auto& cos_theta = get<Tags::detail::CosTheta>(computation_variables);
841  auto& sin_phi = get<Tags::detail::SinPhi>(computation_variables);
842  auto& sin_theta = get<Tags::detail::SinTheta>(computation_variables);
843  trigonometric_functions_on_swsh_collocation(
844  make_not_null(&cos_phi), make_not_null(&cos_theta),
845  make_not_null(&sin_phi), make_not_null(&sin_theta), l_max);
846 
847  // NOTE: to handle the singular values of polar coordinates, the phi
848  // components of all tensors are scaled according to their sin(theta)
849  // prefactors.
850  // so, any down-index component get<2>(A) represents 1/sin(theta) A_\phi,
851  // and any up-index component get<2>(A) represents sin(theta) A^\phi.
852  // This holds for Jacobians, and so direct application of the Jacobians
853  // brings the factors through.
854  auto& cartesian_coords =
855  get<Tags::detail::CartesianCoordinates>(computation_variables);
856  auto& cartesian_to_spherical_jacobian =
857  get<Tags::detail::CartesianToSphericalJacobian>(computation_variables);
858  auto& inverse_cartesian_to_spherical_jacobian =
859  get<Tags::detail::InverseCartesianToSphericalJacobian>(
860  computation_variables);
861  cartesian_to_spherical_coordinates_and_jacobians(
862  make_not_null(&cartesian_coords),
863  make_not_null(&cartesian_to_spherical_jacobian),
864  make_not_null(&inverse_cartesian_to_spherical_jacobian), cos_phi,
865  cos_theta, sin_phi, sin_theta, extraction_radius);
866 
867  auto& spatial_metric =
868  get<gr::Tags::SpatialMetric<3, ::Frame::Inertial, DataVector>>(
869  computation_variables);
871 
872  auto& inverse_spatial_metric =
873  get<gr::Tags::InverseSpatialMetric<3, ::Frame::Inertial, DataVector>>(
874  computation_variables);
875  // Allocation
876  inverse_spatial_metric = determinant_and_inverse(spatial_metric).second;
877 
878  auto& shift = get<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>(
879  computation_variables);
880  gr::shift(make_not_null(&shift), spacetime_metric, inverse_spatial_metric);
881 
882  auto& lapse = get<gr::Tags::Lapse<DataVector>>(computation_variables);
884 
885  auto& dt_spacetime_metric = get<
886  ::Tags::dt<gr::Tags::SpacetimeMetric<3, ::Frame::Inertial, DataVector>>>(
887  computation_variables);
888 
890  make_not_null(&dt_spacetime_metric), lapse, shift, pi, phi);
891 
892  auto& dt_worldtube_normal =
893  get<::Tags::dt<Tags::detail::WorldtubeNormal>>(computation_variables);
894  auto& worldtube_normal =
895  get<Tags::detail::WorldtubeNormal>(computation_variables);
896  worldtube_normal_and_derivatives(
897  make_not_null(&worldtube_normal), make_not_null(&dt_worldtube_normal),
898  cos_phi, cos_theta, spacetime_metric, dt_spacetime_metric, sin_phi,
899  sin_theta, inverse_spatial_metric);
900  auto& spacetime_unit_normal =
901  get<gr::Tags::SpacetimeNormalVector<3, ::Frame::Inertial, DataVector>>(
902  computation_variables);
903  gr::spacetime_normal_vector(make_not_null(&spacetime_unit_normal), lapse,
904  shift);
905  auto& dt_lapse =
906  get<::Tags::dt<gr::Tags::Lapse<DataVector>>>(computation_variables);
908  make_not_null(&dt_lapse), lapse, shift, spacetime_unit_normal, phi, pi);
909  auto& dt_shift =
910  get<::Tags::dt<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>>(
911  computation_variables);
913  shift, inverse_spatial_metric,
914  spacetime_unit_normal, phi, pi);
915 
916  auto& du_null_l = get<::Tags::dt<Tags::detail::NullL>>(computation_variables);
917  auto& null_l = get<Tags::detail::NullL>(computation_variables);
918  null_vector_l_and_derivatives(make_not_null(&du_null_l),
919  make_not_null(&null_l), dt_worldtube_normal,
920  dt_lapse, dt_spacetime_metric, dt_shift, lapse,
921  spacetime_metric, shift, worldtube_normal);
922 
923  // pass to the next step that is common between the 'modal' input and 'GH'
924  // input strategies
925  detail::create_bondi_boundary_data(
926  bondi_boundary_data, make_not_null(&computation_variables),
927  make_not_null(&derivative_buffers), dt_spacetime_metric, phi,
928  spacetime_metric, null_l, du_null_l, cartesian_to_spherical_jacobian,
929  l_max, extraction_radius);
930 }
931 
932 /*!
933  * \brief Process the worldtube data from modal metric components and
934  * derivatives to desired Bondi quantities, placing the result in the passed
935  * \ref DataBoxGroup.
936  *
937  * \details
938  * The mathematics are a bit complicated for all of the coordinate
939  * transformations that are necessary to obtain the Bondi gauge quantities.
940  * For full mathematical details, see the documentation for functions in
941  * `BoundaryData.hpp` and \cite Barkett2019uae \cite Bishop1998uk.
942  *
943  * This function takes as input the full set of ADM metric data and its radial
944  * and time derivatives on a two-dimensional surface of constant \f$r\f$ and
945  * \f$t\f$ in numerical coordinates. This data must be provided as spherical
946  * harmonic coefficients in the libsharp format. This data is provided in nine
947  * `Tensor`s.
948  *
949  * Sufficient tags to provide full worldtube boundary data at a particular
950  * time are set in `bondi_boundary_data`. In particular, the set of tags in
951  * `Tags::characteristic_worldtube_boundary_tags` in the provided \ref
952  * DataBoxGroup are assigned to the worldtube boundary values associated with
953  * the input metric components.
954  *
955  * The majority of the mathematical transformations are implemented as a set of
956  * individual cascaded functions below. The details of the manipulations that
957  * are performed to the input data may be found in the individual functions
958  * themselves, which are called in the following order:
959  * - `trigonometric_functions_on_swsh_collocation()`
960  * - `cartesian_to_spherical_coordinates_and_jacobians()`
961  * - `cartesian_spatial_metric_and_derivatives_from_modes()`
962  * - `cartesian_shift_and_derivatives_from_modes()`
963  * - `cartesian_lapse_and_derivatives_from_modes()`
964  * - `GeneralizedHarmonic::phi()`
965  * - `gr::time_derivative_of_spacetime_metric`
966  * - `gr::spacetime_metric`
967  * - `generalized_harmonic_quantities()`
968  * - `worldtube_normal_and_derivatives()`
969  * - `null_vector_l_and_derivatives()`
970  * - `null_metric_and_derivative()`
971  * - `dlambda_null_metric_and_inverse()`
972  * - `bondi_r()`
973  * - `d_bondi_r()`
974  * - `dyads()`
975  * - `beta_worldtube_data()`
976  * - `bondi_u_worldtube_data()`
977  * - `bondi_w_worldtube_data()`
978  * - `bondi_j_worldtube_data()`
979  * - `dr_bondi_j()`
980  * - `d2lambda_bondi_r()`
981  * - `bondi_q_worldtube_data()`
982  * - `bondi_h_worldtube_data()`
983  * - `du_j_worldtube_data()`
984  */
985 template <typename DataBoxTagList>
987  const gsl::not_null<db::DataBox<DataBoxTagList>*> bondi_boundary_data,
988  const tnsr::ii<ComplexModalVector, 3>& spatial_metric_coefficients,
989  const tnsr::ii<ComplexModalVector, 3>& dt_spatial_metric_coefficients,
990  const tnsr::ii<ComplexModalVector, 3>& dr_spatial_metric_coefficients,
991  const tnsr::I<ComplexModalVector, 3>& shift_coefficients,
992  const tnsr::I<ComplexModalVector, 3>& dt_shift_coefficients,
993  const tnsr::I<ComplexModalVector, 3>& dr_shift_coefficients,
994  const Scalar<ComplexModalVector>& lapse_coefficients,
995  const Scalar<ComplexModalVector>& dt_lapse_coefficients,
996  const Scalar<ComplexModalVector>& dr_lapse_coefficients,
997  const double extraction_radius, const size_t l_max) noexcept {
998  const size_t size = Spectral::Swsh::number_of_swsh_collocation_points(l_max);
999 
1000  // Most allocations required for the full boundary computation are merged into
1001  // a single, large Variables allocation. There remain a handful of cases in
1002  // the computational functions called where an intermediate quantity that is
1003  // not re-used is allocated rather than taking a buffer. These cases are
1004  // marked with code comments 'Allocation'; In the future, if allocations are
1005  // identified as a point to optimize, those buffers may be allocated here and
1006  // passed as function arguments
1007  Variables<tmpl::list<
1008  Tags::detail::CosPhi, Tags::detail::CosTheta, Tags::detail::SinPhi,
1009  Tags::detail::SinTheta, Tags::detail::CartesianCoordinates,
1010  Tags::detail::CartesianToSphericalJacobian,
1011  Tags::detail::InverseCartesianToSphericalJacobian,
1015  tmpl::size_t<3>, ::Frame::Inertial>,
1019  tmpl::size_t<3>, ::Frame::Inertial>,
1023  ::Frame::Inertial>,
1028  Tags::detail::WorldtubeNormal, ::Tags::dt<Tags::detail::WorldtubeNormal>,
1029  Tags::detail::NullL, ::Tags::dt<Tags::detail::NullL>,
1030  // for the detail function called at the end
1034  Tags::detail::AngularDNullL,
1035  Tags::detail::DLambda<
1036  gr::Tags::SpacetimeMetric<3, Frame::RadialNull, DataVector>>,
1037  Tags::detail::DLambda<
1038  gr::Tags::InverseSpacetimeMetric<3, Frame::RadialNull, DataVector>>,
1041  Tags::detail::DLambda<Tags::detail::DLambda<Tags::detail::RealBondiR>>,
1043  tmpl::size_t<2>, Frame::RadialNull>>>
1044  computation_variables{size};
1045 
1046  Variables<
1047  tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
1051  derivative_buffers{size};
1052  auto& cos_phi = get<Tags::detail::CosPhi>(computation_variables);
1053  auto& cos_theta = get<Tags::detail::CosTheta>(computation_variables);
1054  auto& sin_phi = get<Tags::detail::SinPhi>(computation_variables);
1055  auto& sin_theta = get<Tags::detail::SinTheta>(computation_variables);
1056  trigonometric_functions_on_swsh_collocation(
1057  make_not_null(&cos_phi), make_not_null(&cos_theta),
1058  make_not_null(&sin_phi), make_not_null(&sin_theta), l_max);
1059 
1060  // NOTE: to handle the singular values of polar coordinates, the phi
1061  // components of all tensors are scaled according to their sin(theta)
1062  // prefactors.
1063  // so, any down-index component get<2>(A) represents 1/sin(theta) A_\phi,
1064  // and any up-index component get<2>(A) represents sin(theta) A^\phi.
1065  // This holds for Jacobians, and so direct application of the Jacobians
1066  // brings the factors through.
1067  auto& cartesian_coords =
1068  get<Tags::detail::CartesianCoordinates>(computation_variables);
1069  auto& cartesian_to_spherical_jacobian =
1070  get<Tags::detail::CartesianToSphericalJacobian>(computation_variables);
1071  auto& inverse_cartesian_to_spherical_jacobian =
1072  get<Tags::detail::InverseCartesianToSphericalJacobian>(
1073  computation_variables);
1074  cartesian_to_spherical_coordinates_and_jacobians(
1075  make_not_null(&cartesian_coords),
1076  make_not_null(&cartesian_to_spherical_jacobian),
1077  make_not_null(&inverse_cartesian_to_spherical_jacobian), cos_phi,
1078  cos_theta, sin_phi, sin_theta, extraction_radius);
1079 
1080  auto& cartesian_spatial_metric =
1081  get<gr::Tags::SpatialMetric<3, ::Frame::Inertial, DataVector>>(
1082  computation_variables);
1083  auto& inverse_spatial_metric =
1084  get<gr::Tags::InverseSpatialMetric<3, ::Frame::Inertial, DataVector>>(
1085  computation_variables);
1086  auto& d_cartesian_spatial_metric = get<
1087  ::Tags::deriv<gr::Tags::SpatialMetric<3, ::Frame::Inertial, DataVector>,
1088  tmpl::size_t<3>, ::Frame::Inertial>>(computation_variables);
1089  auto& dt_cartesian_spatial_metric = get<
1090  ::Tags::dt<gr::Tags::SpatialMetric<3, ::Frame::Inertial, DataVector>>>(
1091  computation_variables);
1092  auto& interpolation_buffer =
1093  get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
1094  std::integral_constant<int, 0>>>(
1095  derivative_buffers);
1096  Scalar<SpinWeighted<ComplexModalVector, 0>> interpolation_modal_buffer{size};
1097  auto& eth_buffer =
1098  get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
1100  derivative_buffers);
1101  cartesian_spatial_metric_and_derivatives_from_modes(
1102  make_not_null(&cartesian_spatial_metric),
1103  make_not_null(&inverse_spatial_metric),
1104  make_not_null(&d_cartesian_spatial_metric),
1105  make_not_null(&dt_cartesian_spatial_metric),
1106  make_not_null(&interpolation_modal_buffer),
1107  make_not_null(&interpolation_buffer), make_not_null(&eth_buffer),
1108  spatial_metric_coefficients, dr_spatial_metric_coefficients,
1109  dt_spatial_metric_coefficients, inverse_cartesian_to_spherical_jacobian,
1110  l_max);
1111 
1112  auto& cartesian_shift =
1113  get<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>(
1114  computation_variables);
1115  auto& d_cartesian_shift =
1116  get<::Tags::deriv<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>,
1117  tmpl::size_t<3>, ::Frame::Inertial>>(
1118  computation_variables);
1119  auto& dt_cartesian_shift =
1120  get<::Tags::dt<gr::Tags::Shift<3, ::Frame::Inertial, DataVector>>>(
1121  computation_variables);
1122 
1123  cartesian_shift_and_derivatives_from_modes(
1124  make_not_null(&cartesian_shift), make_not_null(&d_cartesian_shift),
1125  make_not_null(&dt_cartesian_shift),
1126  make_not_null(&interpolation_modal_buffer),
1127  make_not_null(&interpolation_buffer), make_not_null(&eth_buffer),
1128  shift_coefficients, dr_shift_coefficients, dt_shift_coefficients,
1129  inverse_cartesian_to_spherical_jacobian, l_max);
1130 
1131  auto& cartesian_lapse =
1132  get<gr::Tags::Lapse<DataVector>>(computation_variables);
1133  auto& d_cartesian_lapse =
1134  get<::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
1135  ::Frame::Inertial>>(computation_variables);
1136  auto& dt_cartesian_lapse =
1137  get<::Tags::dt<gr::Tags::Lapse<DataVector>>>(computation_variables);
1138  cartesian_lapse_and_derivatives_from_modes(
1139  make_not_null(&cartesian_lapse), make_not_null(&d_cartesian_lapse),
1140  make_not_null(&dt_cartesian_lapse),
1141  make_not_null(&interpolation_modal_buffer),
1142  make_not_null(&interpolation_buffer), make_not_null(&eth_buffer),
1143  lapse_coefficients, dr_lapse_coefficients, dt_lapse_coefficients,
1144  inverse_cartesian_to_spherical_jacobian, l_max);
1145 
1146  auto& phi = get<GeneralizedHarmonic::Tags::Phi<3, ::Frame::Inertial>>(
1147  computation_variables);
1148  auto& dt_spacetime_metric = get<
1149  ::Tags::dt<gr::Tags::SpacetimeMetric<3, ::Frame::Inertial, DataVector>>>(
1150  computation_variables);
1151  auto& spacetime_metric =
1152  get<gr::Tags::SpacetimeMetric<3, ::Frame::Inertial, DataVector>>(
1153  computation_variables);
1155  make_not_null(&phi), cartesian_lapse, d_cartesian_lapse, cartesian_shift,
1156  d_cartesian_shift, cartesian_spatial_metric, d_cartesian_spatial_metric);
1158  make_not_null(&dt_spacetime_metric), cartesian_lapse, dt_cartesian_lapse,
1159  cartesian_shift, dt_cartesian_shift, cartesian_spatial_metric,
1160  dt_cartesian_spatial_metric);
1161  gr::spacetime_metric(make_not_null(&spacetime_metric), cartesian_lapse,
1162  cartesian_shift, cartesian_spatial_metric);
1163 
1164  auto& dt_worldtube_normal =
1165  get<::Tags::dt<Tags::detail::WorldtubeNormal>>(computation_variables);
1166  auto& worldtube_normal =
1167  get<Tags::detail::WorldtubeNormal>(computation_variables);
1168  worldtube_normal_and_derivatives(
1169  make_not_null(&worldtube_normal), make_not_null(&dt_worldtube_normal),
1170  cos_phi, cos_theta, spacetime_metric, dt_spacetime_metric, sin_phi,
1171  sin_theta, inverse_spatial_metric);
1172 
1173  auto& du_null_l = get<::Tags::dt<Tags::detail::NullL>>(computation_variables);
1174  auto& null_l = get<Tags::detail::NullL>(computation_variables);
1175  null_vector_l_and_derivatives(
1176  make_not_null(&du_null_l), make_not_null(&null_l), dt_worldtube_normal,
1177  dt_cartesian_lapse, dt_spacetime_metric, dt_cartesian_shift,
1178  cartesian_lapse, spacetime_metric, cartesian_shift, worldtube_normal);
1179 
1180  // pass to the next step that is common between the 'modal' input and 'GH'
1181  // input strategies
1182  detail::create_bondi_boundary_data(
1183  bondi_boundary_data, make_not_null(&computation_variables),
1184  make_not_null(&derivative_buffers), dt_spacetime_metric, phi,
1185  spacetime_metric, null_l, du_null_l, cartesian_to_spherical_jacobian,
1186  l_max, extraction_radius);
1187 }
1188 } // namespace Cce
Auxiliary variable which is analytically the spatial derivative of the spacetime metric.
Definition: Tags.hpp:40
Defines function computing the determinant and inverse of a tensor.
void phi(gsl::not_null< tnsr::iaa< DataType, SpatialDim, Frame > *> phi, const Scalar< DataType > &lapse, const tnsr::i< DataType, SpatialDim, Frame > &deriv_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::iJ< DataType, SpatialDim, Frame > &deriv_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ijj< DataType, SpatialDim, Frame > &deriv_spatial_metric) noexcept
Computes the auxiliary variable used by the generalized harmonic formulation of Einstein&#39;s equations...
Bondi parameter .
Definition: Tags.hpp:77
Definition: Tags.hpp:17
The set of utilities for performing Cauchy characteristic evolution and Cauchy characteristic matchin...
Definition: BoundaryComputeAndSendToEvolution.hpp:28
tnsr::ii< DataType, SpatialDim, Frame > spatial_metric(const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric) noexcept
Compute spatial metric from spacetime metric.
constexpr size_t number_of_swsh_collocation_points(const size_t l_max) noexcept
Convenience function for determining the number of spin-weighted spherical harmonic collocation value...
Definition: SwshCollocation.hpp:25
The derivative with respect to Bondi .
Definition: Tags.hpp:126
Bondi parameter .
Definition: Tags.hpp:34
db::wrap_tags_in< BoundaryPrefix, tmpl::list< Tags::BondiBeta, Tags::BondiU, Tags::Dr< Tags::BondiU >, Tags::BondiQ, Tags::BondiW, Tags::BondiJ, Tags::Dr< Tags::BondiJ >, Tags::BondiH, Tags::Du< Tags::BondiJ >, Tags::BondiR, Tags::Du< Tags::BondiR >, Tags::DuRDividedByR > > characteristic_worldtube_boundary_tags
The collection of tags mutated by create_bondi_boundary_data
Definition: BoundaryData.hpp:545
void time_deriv_of_lapse(gsl::not_null< Scalar< DataType > *> dt_lapse, const Scalar< DataType > &lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::A< DataType, SpatialDim, Frame > &spacetime_unit_normal, const tnsr::iaa< DataType, SpatialDim, Frame > &phi, const tnsr::aa< DataType, SpatialDim, Frame > &pi) noexcept
Computes time derivative of lapse (N) from the generalized harmonic variables, lapse, shift and the spacetime unit normal 1-form.
Define prefixes for DataBox tags.
void create_bondi_boundary_data(const gsl::not_null< db::DataBox< DataBoxTagList > *> bondi_boundary_data, const tnsr::iaa< DataVector, 3 > &phi, const tnsr::aa< DataVector, 3 > &pi, const tnsr::aa< DataVector, 3 > &spacetime_metric, const double extraction_radius, const size_t l_max) noexcept
Process the worldtube data from generalized harmonic quantities to desired Bondi quantities, placing the result in the passed DataBox.
Definition: BoundaryData.hpp:789
The Bondi radius is of the worldtube.
Definition: Tags.hpp:329
tmpl::transform< TagList, tmpl::bind< Wrapper, tmpl::_1, tmpl::pin< Args >... > > wrap_tags_in
Create a new list of Tags by wrapping each tag in TagList using the Wrapper.
Definition: PrefixHelpers.hpp:30
Definition: Tags.hpp:52
Inverse of the spatial metric.
Definition: Tags.hpp:33
Prefix indicating spatial derivatives.
Definition: PartialDerivatives.hpp:54
Defines Functions for calculating spacetime tensors from 3+1 quantities.
void spacetime_metric(gsl::not_null< tnsr::aa< DataType, Dim, Frame > *> spacetime_metric, const Scalar< DataType > &lapse, const tnsr::I< DataType, Dim, Frame > &shift, const tnsr::ii< DataType, Dim, Frame > &spatial_metric) noexcept
Computes the spacetime metric from the spatial metric, lapse, and shift.
Prefix indicating a time derivative.
Definition: Prefixes.hpp:32
Prefix indicating spacetime derivatives.
Definition: PartialDerivatives.hpp:91
Scalar< DataType > lapse(const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric) noexcept
Compute lapse from shift and spacetime metric.
The value , where is Bondi radius of the worldtube.
Definition: Tags.hpp:290
Definition: Determinant.hpp:11
Given a Tag with a type of Tensor<VectorType, ...>, acts as a new version of the tag with type of Ten...
Definition: Tags.hpp:59
Make a spin-weighted type T with spin-weight Spin. Mathematical operators are restricted to addition...
Definition: SpinWeighted.hpp:24
void time_derivative_of_spacetime_metric(gsl::not_null< tnsr::aa< DataType, SpatialDim, Frame > *> dt_spacetime_metric, const Scalar< DataType > &lapse, const Scalar< DataType > &dt_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::I< DataType, SpatialDim, Frame > &dt_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ii< DataType, SpatialDim, Frame > &dt_spatial_metric) noexcept
Computes the time derivative of the spacetime metric from spatial metric, lapse, shift, and their time derivatives.
Definition: Tags.hpp:98
Defines classes and functions used for manipulating DataBox&#39;s.
auto determinant_and_inverse(const Tensor< T, Symm, tmpl::list< Index0, Index1 >> &tensor) noexcept -> std::pair< Scalar< T >, Tensor< T, Symm, tmpl::list< change_index_up_lo< Index1 >, change_index_up_lo< Index0 >>>>
Computes the determinant and inverse of a rank-2 Tensor.
Definition: DeterminantAndInverse.hpp:357
Definition: Tags.hpp:26
The frame for the spherical metric in which the radial coordinate is an affine parameter along outwar...
Definition: BoundaryDataTags.hpp:14
void time_derivative_of_spacetime_metric(gsl::not_null< tnsr::aa< DataType, SpatialDim, Frame > *> dt_spacetime_metric, const Scalar< DataType > &lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::aa< DataType, SpatialDim, Frame > &pi, const tnsr::iaa< DataType, SpatialDim, Frame > &phi) noexcept
Computes the time derivative of the spacetime metric from the generalized harmonic quantities ...
Stores a collection of complex function values.
Definition: ComplexDataVector.hpp:47
Definition: InterpolationTargetWedgeSectionTorus.hpp:25
Definition: Tags.hpp:48
Definition: DataBoxTag.hpp:27
void pi(gsl::not_null< tnsr::aa< DataType, SpatialDim, Frame > *> pi, const Scalar< DataType > &lapse, const Scalar< DataType > &dt_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::I< DataType, SpatialDim, Frame > &dt_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ii< DataType, SpatialDim, Frame > &dt_spatial_metric, const tnsr::iaa< DataType, SpatialDim, Frame > &phi) noexcept
Computes the conjugate momentum of the spacetime metric .
Defines a list of useful type aliases for tensors.
tnsr::I< DataType, SpatialDim, Frame > shift(const tnsr::aa< DataType, SpatialDim, Frame > &spacetime_metric, const tnsr::II< DataType, SpatialDim, Frame > &inverse_spatial_metric) noexcept
Compute shift from spacetime metric and inverse spatial metric.
Stores a collection of function values.
Definition: DataVector.hpp:42
::Tags::dt< BondiJ > BondiH
Bondi parameter .
Definition: Tags.hpp:62
Definition: Tags.hpp:21
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion, but it may be necessary to perform the conversion explicitly when type deduction is desired.
Definition: Gsl.hpp:879
The derivative with respect to Bondi retarded time .
Definition: Tags.hpp:136
Definition: IndexType.hpp:44
void time_deriv_of_shift(gsl::not_null< tnsr::I< DataType, SpatialDim, Frame > *> dt_shift, const Scalar< DataType > &lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::II< DataType, SpatialDim, Frame > &inverse_spatial_metric, const tnsr::A< DataType, SpatialDim, Frame > &spacetime_unit_normal, const tnsr::iaa< DataType, SpatialDim, Frame > &phi, const tnsr::aa< DataType, SpatialDim, Frame > &pi) noexcept
Computes time derivative of the shift vector from the generalized harmonic and geometric variables...
tnsr::A< DataType, SpatialDim, Frame > spacetime_normal_vector(const Scalar< DataType > &lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift) noexcept
Computes spacetime normal vector from lapse and shift.
Defines Functions for calculating spacetime tensors from 3+1 quantities.
A prefix tag representing the boundary data for a quantity on the extraction surface.
Definition: Tags.hpp:204
decltype(auto) constexpr square(const T &x)
Compute the square of x
Definition: ConstantExpressions.hpp:54
Tensor< T, Symmetry<>, index_list<> > Scalar
Scalar type.
Definition: TypeAliases.hpp:21
Bondi parameter .
Definition: Tags.hpp:107
Require a pointer to not be a nullptr
Definition: Gsl.hpp:182