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