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