SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/Cce - BoundaryData.hpp Hit Total Coverage
Commit: 817e13c5144619b701c7cd870655d8dbf94ab8ce Lines: 23 26 88.5 %
Date: 2024-07-19 22:17:05
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : 
       8             : #include "DataStructures/DataBox/DataBox.hpp"
       9             : #include "DataStructures/DataBox/PrefixHelpers.hpp"
      10             : #include "DataStructures/DataBox/Prefixes.hpp"
      11             : #include "DataStructures/DataVector.hpp"
      12             : #include "DataStructures/SpinWeighted.hpp"
      13             : #include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp"
      14             : #include "DataStructures/Tensor/TypeAliases.hpp"
      15             : #include "Evolution/Systems/Cce/BoundaryDataTags.hpp"
      16             : #include "Evolution/Systems/Cce/Tags.hpp"
      17             : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshCollocation.hpp"
      18             : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/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           1 : void trigonometric_functions_on_swsh_collocation(
      46             :     gsl::not_null<Scalar<DataVector>*> cos_phi,
      47             :     gsl::not_null<Scalar<DataVector>*> cos_theta,
      48             :     gsl::not_null<Scalar<DataVector>*> sin_phi,
      49             :     gsl::not_null<Scalar<DataVector>*> sin_theta, size_t l_max);
      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           1 : 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,
      80             :     gsl::not_null<CartesianiSphericalJ*>
      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);
      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           0 : 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,
     102             :     gsl::not_null<Scalar<SpinWeighted<ComplexModalVector, 0>>*>
     103             :         interpolation_modal_buffer,
     104             :     gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
     105             :         interpolation_buffer,
     106             :     gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> eth_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);
     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           1 : 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,
     128             :     gsl::not_null<Scalar<SpinWeighted<ComplexModalVector, 0>>*>
     129             :         interpolation_modal_buffer,
     130             :     gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
     131             :         interpolation_buffer,
     132             :     gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> eth_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);
     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           1 : 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,
     154             :     gsl::not_null<Scalar<SpinWeighted<ComplexModalVector, 0>>*>
     155             :         interpolation_modal_buffer,
     156             :     gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
     157             :         interpolation_buffer,
     158             :     gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> eth_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);
     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           1 : 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);
     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           1 : 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);
     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           1 : 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);
     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           1 : 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);
     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           1 : void bondi_r(gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*> bondi_r,
     279             :              const tnsr::aa<DataVector, 3, Frame::RadialNull>& null_metric);
     280             : 
     281             : /*!
     282             :  * \brief Computes the full 4-dimensional partial of the Bondi radius with
     283             :  * respect to the intermediate null coordinates.
     284             :  *
     285             :  * \details The expression evaluated is obtained from differentiating the
     286             :  * determinant equation for `bondi_r`, from (35) of \cite Barkett2019uae :
     287             :  *
     288             :  * \f[
     289             :  * \partial_\alpha r = \frac{r}{4} \left(g^{A B} \partial_\alpha g_{A B}
     290             :  * - \frac{\partial_\alpha \det q_{A B}}{\det q_{A B}}\right)
     291             :  * \f]
     292             :  *
     293             :  * Note that for the angular derivatives, we just numerically differentiate
     294             :  * using the utilities in `Spectral::Swsh::angular_derivative()`. For the time
     295             :  * and radial derivatives, the second term in the above equation vanishes.
     296             :  */
     297           1 : void d_bondi_r(
     298             :     gsl::not_null<tnsr::a<DataVector, 3, Frame::RadialNull>*> d_bondi_r,
     299             :     const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r,
     300             :     const tnsr::aa<DataVector, 3, Frame::RadialNull>& dlambda_null_metric,
     301             :     const tnsr::aa<DataVector, 3, Frame::RadialNull>& du_null_metric,
     302             :     const tnsr::AA<DataVector, 3, Frame::RadialNull>& inverse_null_metric,
     303             :     size_t l_max);
     304             : 
     305             : /*!
     306             :  * \brief Compute the complex angular dyads used to define the spin-weighted
     307             :  * scalars in the CCE system.
     308             :  *
     309             :  * \details We use the typically chosen angular dyads in CCE
     310             :  * \cite Barkett2019uae \cite Bishop1997ik :
     311             :  *
     312             :  * \f{align*}{
     313             :  * q_A &= \{-1, -i \sin(\theta)\}\\
     314             :  * q^A &= \left\{-1, -i \frac{1}{\sin \theta}\right\}
     315             :  * \f}
     316             :  *
     317             :  * However, to maintain regularity and for compatibility with the more regular
     318             :  * Jacobians from `Cce::cartesian_to_spherical_coordinates_and_jacobians()`, in
     319             :  * the code we omit the factors of \f$\sin \theta\f$ from the above equations.
     320             :  */
     321           1 : void dyads(
     322             :     gsl::not_null<tnsr::i<ComplexDataVector, 2, Frame::RadialNull>*> down_dyad,
     323             :     gsl::not_null<tnsr::I<ComplexDataVector, 2, Frame::RadialNull>*> up_dyad);
     324             : 
     325             : /*!
     326             :  * \brief Compute the \f$\beta\f$ (lapse) function for the CCE Bondi-like
     327             :  * metric.
     328             :  *
     329             :  * \details The Bondi-like metric has \f$g^{u r} = - e^{2 \beta}\f$, and the
     330             :  * value of \f$\beta\f$ is obtained from the intermediate null metric by (see
     331             :  * equation (51) of \cite Barkett2019uae) using:
     332             :  *
     333             :  * \f[
     334             :  * \beta = -\frac{1}{2} \ln \partial_{\lambda} r
     335             :  * \f]
     336             :  */
     337           1 : void beta_worldtube_data(
     338             :     gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*> beta,
     339             :     const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r);
     340             : 
     341             : /*!
     342             :  * \brief Compute the \f$U\f$ (shift) function for the CCE Bondi-like metric.
     343             :  *
     344             :  * \details The Bondi-like metric has \f$g^{r A} = -e^{-2 \beta} U^A\f$, and the
     345             :  * spin-weighted vector \f$U = U^A q_A\f$. The value of \f$U^A\f$ can be
     346             :  * computed from the intermediate null metric quantities (see equation (54) of
     347             :  * \cite Barkett2019uae) using:
     348             :  *
     349             :  * \f[
     350             :  * U = -(\partial_\lambda r g^{\lambda A} + \partial_B r g^{A B}) q_A
     351             :  * / \partial_\lambda r \f]
     352             :  *
     353             :  */
     354           1 : void bondi_u_worldtube_data(
     355             :     gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> bondi_u,
     356             :     const tnsr::i<ComplexDataVector, 2, Frame::RadialNull>& dyad,
     357             :     const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
     358             :     const tnsr::AA<DataVector, 3, Frame::RadialNull>& inverse_null_metric);
     359             : 
     360             : /*!
     361             :  * \brief Compute the \f$W\f$ (mass aspect) function for the CCE Bondi-like
     362             :  * metric.
     363             :  *
     364             :  * \details The Bondi-like metric has \f$g^{rr} = e^{-2 \beta}(1 + r W)\f$. The
     365             :  * value of \f$W\f$ can be computed from the null metric quantities (see
     366             :  * equation (55) of \cite Barkett2019uae) using:
     367             :  *
     368             :  * \f[
     369             :  * W = \frac{1}{r} \left(-1
     370             :  * + \frac{g^{\lambda \lambda} (\partial_\lambda r)^2
     371             :  * + 2 \partial_\lambda r \left(\partial_A r g^{\lambda A}
     372             :  * - \partial_u r\right) + \partial_A r \partial_B r g^{A B}}
     373             :  * {\partial_\lambda r}\right) \f]
     374             :  */
     375           1 : void bondi_w_worldtube_data(
     376             :     gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*> bondi_w,
     377             :     const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
     378             :     const tnsr::AA<DataVector, 3, Frame::RadialNull>& inverse_null_metric,
     379             :     const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r);
     380             : 
     381             : /*!
     382             :  * \brief Compute the \f$J\f$ (intuitively similar to the transverse-traceless
     383             :  * part of the angular metric) function for the CCE Bondi-like metric.
     384             :  *
     385             :  * \details The Bondi-like metric has \f$J = \frac{1}{2 r^2} q^A q^B g_{A B}\f$.
     386             :  * This expression holds both for the right-hand side in the Bondi coordinates
     387             :  * and for the right-hand side in the intermediate null coordinates (see
     388             :  * equation (45) of \cite Barkett2019uae).
     389             :  */
     390           1 : void bondi_j_worldtube_data(
     391             :     gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> bondi_j,
     392             :     const tnsr::aa<DataVector, 3, Frame::RadialNull>& null_metric,
     393             :     const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r,
     394             :     const tnsr::I<ComplexDataVector, 2, Frame::RadialNull>& dyad);
     395             : 
     396             : /*!
     397             :  * \brief Compute the radial derivative of the angular metric spin-weighted
     398             :  * scalar \f$\partial_r J\f$ in the CCE Bondi-like metric.
     399             :  *
     400             :  * \details The radial derivative of the angular spin-weighted scalar \f$J\f$
     401             :  * can be computed from the null metric components by (c.f. equation (47) of
     402             :  * \cite Barkett2019uae):
     403             :  *
     404             :  * \f[
     405             :  * \partial_r J = \frac{\partial_\lambda J}{\partial_\lambda r} =
     406             :  *  \frac{q^A q^B \partial_\lambda g_{A B} / (2 r^2)
     407             :  * - 2 \partial_\lambda r J / r}{\partial_\lambda r}
     408             :  * \f]
     409             :  */
     410           1 : void dr_bondi_j(
     411             :     gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> dr_bondi_j,
     412             :     gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
     413             :         denominator_buffer,
     414             :     const tnsr::aa<DataVector, 3, Frame::RadialNull>& dlambda_null_metric,
     415             :     const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
     416             :     const Scalar<SpinWeighted<ComplexDataVector, 2>>& bondi_j,
     417             :     const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r,
     418             :     const tnsr::I<ComplexDataVector, 2, Frame::RadialNull>& dyad);
     419             : 
     420             : /*!
     421             :  * \brief Compute the second derivative of the Bondi radius with respect to the
     422             :  * intermediate null coordinate radius \f$\partial_\lambda^2 r\f$.
     423             :  *
     424             :  * \details To determine this second derivative quantity without resorting to
     425             :  * depending on second-derivative metric inputs, we need to take advantage of
     426             :  * one of the Einstein field equations. Combining equations (53) and (52) of
     427             :  * \cite Barkett2019uae, we have:
     428             :  *
     429             :  * \f[
     430             :  * \partial_\lambda^2 r = \frac{-r}{4} \left(
     431             :  * \partial_\lambda J \partial_\lambda \bar J - (\partial_\lambda K)^2\right)
     432             :  * \f],
     433             :  *
     434             :  * where the first derivative of \f$K\f$ can be obtained from \f$K = \sqrt{1 + J
     435             :  * \bar J}\f$ and the first derivative of \f$J\f$ can be obtained from (47) of
     436             :  * \cite Barkett2019uae
     437             :  */
     438           1 : void d2lambda_bondi_r(
     439             :     gsl::not_null<Scalar<DataVector>*> d2lambda_bondi_r,
     440             :     const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
     441             :     const Scalar<SpinWeighted<ComplexDataVector, 2>>& dr_bondi_j,
     442             :     const Scalar<SpinWeighted<ComplexDataVector, 2>>& bondi_j,
     443             :     const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r);
     444             : 
     445             : /*!
     446             :  * \brief Compute the Bondi metric contribution \f$Q\f$ (radial derivative of
     447             :  * shift).
     448             :  *
     449             :  * \details The definition of \f$Q\f$ in terms of the Bondi metric components is
     450             :  *
     451             :  * \f[
     452             :  *  Q = q^A e^{-2 \beta} g_{A B} \partial_r U^B.
     453             :  * \f]
     454             :  *
     455             :  * $Q$ can be derived from the intermediate null metric quantities via (see
     456             :  * equations (56) and (57) of \cite Barkett2019uae)
     457             :  *
     458             :  * \f[
     459             :  * \partial_\lambda U = - \left(\partial_\lambda g^{\lambda A}
     460             :  * + \frac{\partial_A \partial_\lambda r}{\partial_\lambda r} g^{A B}
     461             :  * + \frac{\partial_B r}{\partial_\lambda r} \partial_\lambda g^{A B}\right) q_A
     462             :  * + 2 \partial_\lambda \beta (U + g^{\lambda A} q_A)
     463             :  * \f]
     464             :  *
     465             :  * and
     466             :  *
     467             :  * \f[
     468             :  * Q = r^2 (J \partial_\lambda \bar U + K \partial_\lambda U)
     469             :  * \f]
     470             :  *
     471             :  * also provided is \f$\partial_r U\f$, which is separately useful to cache for
     472             :  * other intermediate steps in the CCE computation.
     473             :  */
     474           1 : void bondi_q_worldtube_data(
     475             :     gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> bondi_q,
     476             :     gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> dr_bondi_u,
     477             :     const Scalar<DataVector>& d2lambda_r,
     478             :     const tnsr::AA<DataVector, 3, Frame::RadialNull>&
     479             :         dlambda_inverse_null_metric,
     480             :     const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
     481             :     const tnsr::i<ComplexDataVector, 2, Frame::RadialNull>& dyad,
     482             :     const tnsr::i<DataVector, 2, Frame::RadialNull>& angular_d_dlambda_r,
     483             :     const tnsr::AA<DataVector, 3, Frame::RadialNull>& inverse_null_metric,
     484             :     const Scalar<SpinWeighted<ComplexDataVector, 2>>& bondi_j,
     485             :     const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r,
     486             :     const Scalar<SpinWeighted<ComplexDataVector, 1>>& bondi_u);
     487             : 
     488             : /*!
     489             :  * \brief Compute the Bondi metric contribution \f$(\partial_u J)_{y} \equiv
     490             :  * H\f$ (the retarded time derivative evaluated at fixed $y$ coordinate) on the
     491             :  * worldtube boundary.
     492             :  *
     493             :  * \details The numerical time derivative (along the worldtube, rather than
     494             :  * along the surface of constant Bondi \f$r\f$) is computed by (see equation
     495             :  * (48) of \cite Barkett2019uae)
     496             :  *
     497             :  * \f[
     498             :  * (\partial_u J)_y = \frac{1}{2 r^2} q^A q^B \partial_u g_{A B}
     499             :  * - \frac{2 \partial_u r}{r} J
     500             :  * \f]
     501             :  *
     502             :  * \note There is the regrettable notation difference with the primary reference
     503             :  * for these formulas \cite Barkett2019uae in that we denote with \f$H\f$ the
     504             :  * time derivative at constant numerical radius, where \cite Barkett2019uae uses
     505             :  * \f$H\f$ to denote the time derivative at constant Bondi radius.
     506             :  */
     507           1 : void bondi_h_worldtube_data(
     508             :     gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> bondi_h,
     509             :     const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
     510             :     const Scalar<SpinWeighted<ComplexDataVector, 2>>& bondi_j,
     511             :     const tnsr::aa<DataVector, 3, Frame::RadialNull>& du_null_metric,
     512             :     const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r,
     513             :     const tnsr::I<ComplexDataVector, 2, Frame::RadialNull>& dyad);
     514             : 
     515             : /*!
     516             :  * \brief Compute the Bondi metric contribution \f$(\partial_u J)_r\f$ (the
     517             :  * retarded time derivative at fixed coordinate $r$) on the worldtube boundary.
     518             :  *
     519             :  * \details The numerical time derivative (along the surface of constant r, not
     520             :  * along the worldtube) is computed by (see equation (50) of
     521             :  * \cite Barkett2019uae)
     522             :  *
     523             :  * \f[
     524             :  * \partial_u J = \frac{1}{2 r^2} q^A q^B \left(\partial_u g_{A B}
     525             :  * - \frac{ \partial_u r}{ \partial_\lambda r} \partial_\lambda g_{A B}\right)
     526             :  * \f]
     527             :  *
     528             :  * \note There is the regrettable notation difference with the primary reference
     529             :  * for these formulas \cite Barkett2019uae in that we denote with \f$H\f$ the
     530             :  * time derivative at constant numerical radius, where \cite Barkett2019uae uses
     531             :  * \f$H\f$ to denote the time derivative at constant Bondi radius.
     532             :  */
     533           1 : void du_j_worldtube_data(
     534             :     gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> du_bondi_j,
     535             :     const tnsr::a<DataVector, 3, Frame::RadialNull>& d_bondi_r,
     536             :     const Scalar<SpinWeighted<ComplexDataVector, 2>>& bondi_j,
     537             :     const tnsr::aa<DataVector, 3, Frame::RadialNull>& du_null_metric,
     538             :     const tnsr::aa<DataVector, 3, Frame::RadialNull>& dlambda_null_metric,
     539             :     const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r,
     540             :     const tnsr::I<ComplexDataVector, 2, Frame::RadialNull>& dyad);
     541             : 
     542             : namespace Tags {
     543             : /// The collection of tags mutated by `create_bondi_boundary_data`
     544             : template <template <typename> class BoundaryPrefix>
     545           1 : using characteristic_worldtube_boundary_tags = db::wrap_tags_in<
     546             :     BoundaryPrefix,
     547             :     tmpl::list<Tags::BondiBeta, Tags::BondiU, Tags::Dr<Tags::BondiU>,
     548             :                Tags::BondiQ, Tags::BondiW, Tags::BondiJ, Tags::Dr<Tags::BondiJ>,
     549             :                Tags::BondiH, Tags::Du<Tags::BondiJ>, Tags::BondiR,
     550             :                Tags::Du<Tags::BondiR>, Tags::DuRDividedByR>>;
     551             : 
     552           0 : using klein_gordon_worldtube_boundary_tags =
     553             :     tmpl::list<Tags::BoundaryValue<Tags::KleinGordonPsi>,
     554             :                Tags::BoundaryValue<Tags::KleinGordonPi>>;
     555             : }  // namespace Tags
     556             : 
     557             : namespace detail {
     558             : // the common step between the modal input and the Generalized harmonic input
     559             : // that performs the final gauge processing to Bondi scalars and places them in
     560             : // the Variables.
     561             : template <typename BoundaryTagList, typename BufferTagList,
     562             :           typename ComplexBufferTagList>
     563             : void create_bondi_boundary_data(
     564             :     const gsl::not_null<Variables<BoundaryTagList>*> bondi_boundary_data,
     565             :     const gsl::not_null<Variables<BufferTagList>*> computation_variables,
     566             :     const gsl::not_null<Variables<ComplexBufferTagList>*> derivative_buffers,
     567             :     const tnsr::aa<DataVector, 3>& dt_spacetime_metric,
     568             :     const tnsr::iaa<DataVector, 3>& phi,
     569             :     const tnsr::aa<DataVector, 3>& spacetime_metric,
     570             :     const tnsr::A<DataVector, 3>& null_l,
     571             :     const tnsr::A<DataVector, 3>& du_null_l,
     572             :     const SphericaliCartesianJ& cartesian_to_spherical_jacobian,
     573             :     const size_t l_max, const double extraction_radius) {
     574             :   const size_t size = Spectral::Swsh::number_of_swsh_collocation_points(l_max);
     575             : 
     576             :   // unfortunately, because the dyads are not themselves spin-weighted, they
     577             :   // need a separate Variables
     578             :   Variables<tmpl::list<Tags::detail::DownDyad, Tags::detail::UpDyad>>
     579             :       dyad_variables{size};
     580             : 
     581             :   auto& null_metric =
     582             :       get<gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>>(
     583             :           *computation_variables);
     584             :   auto& du_null_metric = get<
     585             :       ::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>>>(
     586             :       *computation_variables);
     587             :   null_metric_and_derivative(
     588             :       make_not_null(&du_null_metric), make_not_null(&null_metric),
     589             :       cartesian_to_spherical_jacobian, dt_spacetime_metric, spacetime_metric);
     590             : 
     591             :   auto& inverse_null_metric =
     592             :       get<gr::Tags::InverseSpacetimeMetric<DataVector, 3, Frame::RadialNull>>(
     593             :           *computation_variables);
     594             : 
     595             :   // the below scaling process is used to reduce accumulation of numerical
     596             :   // error in the determinant evaluation
     597             : 
     598             :   // buffer reuse because the scaled null metric is only needed until the
     599             :   // `determinant_and_inverse` call
     600             :   auto& scaled_null_metric =
     601             :       get<gr::Tags::InverseSpacetimeMetric<DataVector, 3, Frame::RadialNull>>(
     602             :           *computation_variables);
     603             :   for (size_t i = 0; i < 4; ++i) {
     604             :     for (size_t j = i; j < 4; ++j) {
     605             :       if (i > 1 and j > 1) {
     606             :         scaled_null_metric.get(i, j) =
     607             :             null_metric.get(i, j) / square(extraction_radius);
     608             :       } else if (i > 1 or j > 1) {
     609             :         scaled_null_metric.get(i, j) =
     610             :             null_metric.get(i, j) / extraction_radius;
     611             :       } else {
     612             :         scaled_null_metric.get(i, j) = null_metric.get(i, j);
     613             :       }
     614             :     }
     615             :   }
     616             :   // Allocation
     617             :   const auto scaled_inverse_null_metric =
     618             :       determinant_and_inverse(scaled_null_metric).second;
     619             :   for (size_t i = 0; i < 4; ++i) {
     620             :     for (size_t j = i; j < 4; ++j) {
     621             :       if (i > 1 and j > 1) {
     622             :         inverse_null_metric.get(i, j) =
     623             :             scaled_inverse_null_metric.get(i, j) / square(extraction_radius);
     624             :       } else if (i > 1 or j > 1) {
     625             :         inverse_null_metric.get(i, j) =
     626             :             scaled_inverse_null_metric.get(i, j) / extraction_radius;
     627             :       } else {
     628             :         inverse_null_metric.get(i, j) = scaled_inverse_null_metric.get(i, j);
     629             :       }
     630             :     }
     631             :   }
     632             : 
     633             :   auto& angular_d_null_l =
     634             :       get<Tags::detail::AngularDNullL>(*computation_variables);
     635             :   auto& buffer_for_derivatives =
     636             :       get(get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
     637             :                                    std::integral_constant<int, 0>>>(
     638             :           *derivative_buffers));
     639             :   auto& eth_buffer =
     640             :       get(get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
     641             :                                    std::integral_constant<int, 1>>>(
     642             :           *derivative_buffers));
     643             :   for (size_t a = 0; a < 4; ++a) {
     644             :     buffer_for_derivatives.data() =
     645             :         std::complex<double>(1.0, 0.0) * null_l.get(a);
     646             :     Spectral::Swsh::angular_derivatives<tmpl::list<Spectral::Swsh::Tags::Eth>>(
     647             :         l_max, 1, make_not_null(&eth_buffer), buffer_for_derivatives);
     648             :     angular_d_null_l.get(0, a) = -real(eth_buffer.data());
     649             :     angular_d_null_l.get(1, a) = -imag(eth_buffer.data());
     650             :   }
     651             : 
     652             :   auto& dlambda_null_metric = get<Tags::detail::DLambda<
     653             :       gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>>>(
     654             :       *computation_variables);
     655             :   auto& dlambda_inverse_null_metric = get<Tags::detail::DLambda<
     656             :       gr::Tags::InverseSpacetimeMetric<DataVector, 3, Frame::RadialNull>>>(
     657             :       *computation_variables);
     658             :   dlambda_null_metric_and_inverse(
     659             :       make_not_null(&dlambda_null_metric),
     660             :       make_not_null(&dlambda_inverse_null_metric), angular_d_null_l,
     661             :       cartesian_to_spherical_jacobian, phi, dt_spacetime_metric, du_null_l,
     662             :       inverse_null_metric, null_l, spacetime_metric);
     663             : 
     664             :   auto& r = get<Tags::BoundaryValue<Tags::BondiR>>(*bondi_boundary_data);
     665             :   bondi_r(make_not_null(&r), 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             :   get(get<Tags::BoundaryValue<Tags::DuRDividedByR>>(*bondi_boundary_data))
     673             :       .data() =
     674             :       std::complex<double>(1.0, 0.0) * get<0>(d_r) / get(r).data();
     675             :   get(get<Tags::BoundaryValue<Tags::Du<Tags::BondiR>>>(*bondi_boundary_data))
     676             :       .data() = std::complex<double>(1.0, 0.0) * get<0>(d_r);
     677             : 
     678             :   auto& down_dyad = get<Tags::detail::DownDyad>(dyad_variables);
     679             :   auto& up_dyad = get<Tags::detail::UpDyad>(dyad_variables);
     680             :   dyads(make_not_null(&down_dyad), make_not_null(&up_dyad));
     681             : 
     682             :   beta_worldtube_data(make_not_null(&get<Tags::BoundaryValue<Tags::BondiBeta>>(
     683             :                           *bondi_boundary_data)),
     684             :                       d_r);
     685             : 
     686             :   auto& bondi_u = get<Tags::BoundaryValue<Tags::BondiU>>(*bondi_boundary_data);
     687             :   bondi_u_worldtube_data(make_not_null(&bondi_u), down_dyad, d_r,
     688             :                          inverse_null_metric);
     689             : 
     690             :   bondi_w_worldtube_data(make_not_null(&get<Tags::BoundaryValue<Tags::BondiW>>(
     691             :                              *bondi_boundary_data)),
     692             :                          d_r, inverse_null_metric, r);
     693             : 
     694             :   auto& bondi_j =
     695             :       get<Tags::BoundaryValue<Tags::BondiJ>>(*bondi_boundary_data);
     696             :   bondi_j_worldtube_data(make_not_null(&bondi_j), null_metric, r, up_dyad);
     697             : 
     698             :   auto& dr_j =
     699             :       get<Tags::BoundaryValue<Tags::Dr<Tags::BondiJ>>>(*bondi_boundary_data);
     700             :   auto& denominator_buffer =
     701             :       get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
     702             :                                std::integral_constant<int, 0>>>(
     703             :           *derivative_buffers);
     704             :   dr_bondi_j(make_not_null(&get<Tags::BoundaryValue<Tags::Dr<Tags::BondiJ>>>(
     705             :                  *bondi_boundary_data)),
     706             :              make_not_null(&denominator_buffer), dlambda_null_metric, d_r,
     707             :              bondi_j, r, up_dyad);
     708             : 
     709             :   auto& d2lambda_r = get<
     710             :       Tags::detail::DLambda<Tags::detail::DLambda<Tags::detail::RealBondiR>>>(
     711             :       *computation_variables);
     712             :   d2lambda_bondi_r(make_not_null(&d2lambda_r), d_r, dr_j, bondi_j, r);
     713             : 
     714             :   auto& angular_d_dlambda_r =
     715             :       get<::Tags::deriv<Tags::detail::DLambda<Tags::detail::RealBondiR>,
     716             :                         tmpl::size_t<2>, Frame::RadialNull>>(
     717             :           *computation_variables);
     718             :   buffer_for_derivatives.data() = std::complex<double>(1.0, 0.0) * get<1>(d_r);
     719             :   Spectral::Swsh::angular_derivatives<tmpl::list<Spectral::Swsh::Tags::Eth>>(
     720             :       l_max, 1, make_not_null(&eth_buffer), buffer_for_derivatives);
     721             :   angular_d_dlambda_r.get(0) = -real(eth_buffer.data());
     722             :   angular_d_dlambda_r.get(1) = -imag(eth_buffer.data());
     723             : 
     724             :   bondi_q_worldtube_data(
     725             :       make_not_null(
     726             :           &get<Tags::BoundaryValue<Tags::BondiQ>>(*bondi_boundary_data)),
     727             :       make_not_null(&get<Tags::BoundaryValue<Tags::Dr<Tags::BondiU>>>(
     728             :           *bondi_boundary_data)),
     729             :       d2lambda_r, dlambda_inverse_null_metric, d_r, down_dyad,
     730             :       angular_d_dlambda_r, inverse_null_metric, bondi_j, r, bondi_u);
     731             : 
     732             :   bondi_h_worldtube_data(make_not_null(&get<Tags::BoundaryValue<Tags::BondiH>>(
     733             :                              *bondi_boundary_data)),
     734             :                          d_r, bondi_j, du_null_metric, r, up_dyad);
     735             : 
     736             :   du_j_worldtube_data(
     737             :       make_not_null(&get<Tags::BoundaryValue<Tags::Du<Tags::BondiJ>>>(
     738             :           *bondi_boundary_data)),
     739             :       d_r, bondi_j, du_null_metric, dlambda_null_metric, r, up_dyad);
     740             : }
     741             : }  // namespace detail
     742             : 
     743             : /*!
     744             :  * \brief Process the worldtube data from generalized harmonic quantities
     745             :  *  to desired Bondi quantities, placing the result in the passed
     746             :  * `Variables`.
     747             :  *
     748             :  * \details
     749             :  * The mathematics are a bit complicated for all of the coordinate
     750             :  * transformations that are necessary to obtain the Bondi gauge quantities.
     751             :  * For full mathematical details, see the documentation for functions in
     752             :  * `BoundaryData.hpp` and \cite Barkett2019uae \cite Bishop1998uk.
     753             :  *
     754             :  * This function takes as input the full set of Generalized harmonic metric data
     755             :  * on a two-dimensional surface of constant \f$r\f$ and \f$t\f$ in numerical
     756             :  * coordinates.
     757             :  *
     758             :  * Sufficient tags to provide full worldtube boundary data at a particular
     759             :  * time are set in `bondi_boundary_data`. In particular, the set of tags in
     760             :  * `Tags::characteristic_worldtube_boundary_tags` in the provided `Variables`
     761             :  * are assigned to the worldtube boundary values associated with the input
     762             :  * metric components.
     763             :  *
     764             :  * The majority of the mathematical transformations are implemented as a set of
     765             :  * individual cascaded functions below. The details of the manipulations that
     766             :  * are performed to the input data may be found in the individual functions
     767             :  * themselves, which are called in the following order:
     768             :  * - `trigonometric_functions_on_swsh_collocation()`
     769             :  * - `gr::shift()`
     770             :  * - `gr::lapse()`
     771             :  * - `worldtube_normal_and_derivatives()`
     772             :  * - `gr::spacetime_normal_vector()`
     773             :  * - `gh::time_deriv_of_lapse()`
     774             :  * - `gh::time_deriv_of_shift()`
     775             :  * - `null_vector_l_and_derivatives()`
     776             :  * - `cartesian_to_spherical_coordinates_and_jacobians()`
     777             :  * - `null_metric_and_derivative()`
     778             :  * - `dlambda_null_metric_and_inverse()`
     779             :  * - `bondi_r()`
     780             :  * - `d_bondi_r()`
     781             :  * - `dyads()`
     782             :  * - `beta_worldtube_data()`
     783             :  * - `bondi_u_worldtube_data()`
     784             :  * - `bondi_w_worldtube_data()`
     785             :  * - `bondi_j_worldtube_data()`
     786             :  * - `dr_bondi_j()`
     787             :  * - `d2lambda_bondi_r()`
     788             :  * - `bondi_q_worldtube_data()`
     789             :  * - `bondi_h_worldtube_data()`
     790             :  * - `du_j_worldtube_data()`
     791             :  */
     792             : template <typename BoundaryTagList>
     793           1 : void create_bondi_boundary_data(
     794             :     const gsl::not_null<Variables<BoundaryTagList>*> bondi_boundary_data,
     795             :     const tnsr::iaa<DataVector, 3>& phi, const tnsr::aa<DataVector, 3>& pi,
     796             :     const tnsr::aa<DataVector, 3>& spacetime_metric,
     797             :     const double extraction_radius, const size_t l_max) {
     798             :   const size_t size = Spectral::Swsh::number_of_swsh_collocation_points(l_max);
     799             :   // Most allocations required for the full boundary computation are merged into
     800             :   // a single, large Variables allocation. There remain a handful of cases in
     801             :   // the computational functions called where an intermediate quantity that is
     802             :   // not re-used is allocated rather than taking a buffer. These cases are
     803             :   // marked with code comments 'Allocation'; In the future, if allocations are
     804             :   // identified as a point to optimize, those buffers may be allocated here and
     805             :   // passed as function arguments
     806             :   Variables<tmpl::list<
     807             :       Tags::detail::CosPhi, Tags::detail::CosTheta, Tags::detail::SinPhi,
     808             :       Tags::detail::SinTheta, Tags::detail::CartesianCoordinates,
     809             :       Tags::detail::CartesianToSphericalJacobian,
     810             :       Tags::detail::InverseCartesianToSphericalJacobian,
     811             :       gr::Tags::SpatialMetric<DataVector, 3>,
     812             :       gr::Tags::InverseSpatialMetric<DataVector, 3>,
     813             :       gr::Tags::Shift<DataVector, 3>,
     814             :       ::Tags::dt<gr::Tags::Shift<DataVector, 3>>, gr::Tags::Lapse<DataVector>,
     815             :       ::Tags::dt<gr::Tags::Lapse<DataVector>>,
     816             :       ::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3>>,
     817             :       Tags::detail::WorldtubeNormal, ::Tags::dt<Tags::detail::WorldtubeNormal>,
     818             :       gr::Tags::SpacetimeNormalVector<DataVector, 3>, Tags::detail::NullL,
     819             :       ::Tags::dt<Tags::detail::NullL>,
     820             :       // for the detail function called at the end
     821             :       gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>,
     822             :       ::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>>,
     823             :       gr::Tags::InverseSpacetimeMetric<DataVector, 3, Frame::RadialNull>,
     824             :       Tags::detail::AngularDNullL,
     825             :       Tags::detail::DLambda<
     826             :           gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>>,
     827             :       Tags::detail::DLambda<
     828             :           gr::Tags::InverseSpacetimeMetric<DataVector, 3, Frame::RadialNull>>,
     829             :       ::Tags::spacetime_deriv<Tags::detail::RealBondiR, tmpl::size_t<3>,
     830             :                               Frame::RadialNull>,
     831             :       Tags::detail::DLambda<Tags::detail::DLambda<Tags::detail::RealBondiR>>,
     832             :       ::Tags::deriv<Tags::detail::DLambda<Tags::detail::RealBondiR>,
     833             :                     tmpl::size_t<2>, Frame::RadialNull>>>
     834             :       computation_variables{size};
     835             : 
     836             :   Variables<
     837             :       tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
     838             :                                       std::integral_constant<int, 0>>,
     839             :                  ::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
     840             :                                       std::integral_constant<int, 1>>>>
     841             :       derivative_buffers{size};
     842             : 
     843             :   auto& cos_phi = get<Tags::detail::CosPhi>(computation_variables);
     844             :   auto& cos_theta = get<Tags::detail::CosTheta>(computation_variables);
     845             :   auto& sin_phi = get<Tags::detail::SinPhi>(computation_variables);
     846             :   auto& sin_theta = get<Tags::detail::SinTheta>(computation_variables);
     847             :   trigonometric_functions_on_swsh_collocation(
     848             :       make_not_null(&cos_phi), make_not_null(&cos_theta),
     849             :       make_not_null(&sin_phi), make_not_null(&sin_theta), l_max);
     850             : 
     851             :   // NOTE: to handle the singular values of polar coordinates, the phi
     852             :   // components of all tensors are scaled according to their sin(theta)
     853             :   // prefactors.
     854             :   // so, any down-index component get<2>(A) represents 1/sin(theta) A_\phi,
     855             :   // and any up-index component get<2>(A) represents sin(theta) A^\phi.
     856             :   // This holds for Jacobians, and so direct application of the Jacobians
     857             :   // brings the factors through.
     858             :   auto& cartesian_coords =
     859             :       get<Tags::detail::CartesianCoordinates>(computation_variables);
     860             :   auto& cartesian_to_spherical_jacobian =
     861             :       get<Tags::detail::CartesianToSphericalJacobian>(computation_variables);
     862             :   auto& inverse_cartesian_to_spherical_jacobian =
     863             :       get<Tags::detail::InverseCartesianToSphericalJacobian>(
     864             :           computation_variables);
     865             :   cartesian_to_spherical_coordinates_and_jacobians(
     866             :       make_not_null(&cartesian_coords),
     867             :       make_not_null(&cartesian_to_spherical_jacobian),
     868             :       make_not_null(&inverse_cartesian_to_spherical_jacobian), cos_phi,
     869             :       cos_theta, sin_phi, sin_theta, extraction_radius);
     870             : 
     871             :   auto& spatial_metric =
     872             :       get<gr::Tags::SpatialMetric<DataVector, 3>>(computation_variables);
     873             :   gr::spatial_metric(make_not_null(&spatial_metric), spacetime_metric);
     874             : 
     875             :   auto& inverse_spatial_metric =
     876             :       get<gr::Tags::InverseSpatialMetric<DataVector, 3>>(computation_variables);
     877             :   // Allocation
     878             :   inverse_spatial_metric = determinant_and_inverse(spatial_metric).second;
     879             : 
     880             :   auto& shift = get<gr::Tags::Shift<DataVector, 3>>(computation_variables);
     881             :   gr::shift(make_not_null(&shift), spacetime_metric, inverse_spatial_metric);
     882             : 
     883             :   auto& lapse = get<gr::Tags::Lapse<DataVector>>(computation_variables);
     884             :   gr::lapse(make_not_null(&lapse), shift, spacetime_metric);
     885             : 
     886             :   auto& dt_spacetime_metric =
     887             :       get<::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3>>>(
     888             :           computation_variables);
     889             : 
     890             :   gh::time_derivative_of_spacetime_metric(make_not_null(&dt_spacetime_metric),
     891             :                                           lapse, shift, pi, phi);
     892             : 
     893             :   auto& dt_worldtube_normal =
     894             :       get<::Tags::dt<Tags::detail::WorldtubeNormal>>(computation_variables);
     895             :   auto& worldtube_normal =
     896             :       get<Tags::detail::WorldtubeNormal>(computation_variables);
     897             :   worldtube_normal_and_derivatives(
     898             :       make_not_null(&worldtube_normal), make_not_null(&dt_worldtube_normal),
     899             :       cos_phi, cos_theta, spacetime_metric, dt_spacetime_metric, sin_phi,
     900             :       sin_theta, inverse_spatial_metric);
     901             :   auto& spacetime_unit_normal =
     902             :       get<gr::Tags::SpacetimeNormalVector<DataVector, 3>>(
     903             :           computation_variables);
     904             :   gr::spacetime_normal_vector(make_not_null(&spacetime_unit_normal), lapse,
     905             :                               shift);
     906             :   auto& dt_lapse =
     907             :       get<::Tags::dt<gr::Tags::Lapse<DataVector>>>(computation_variables);
     908             :   gh::time_deriv_of_lapse(make_not_null(&dt_lapse), lapse, shift,
     909             :                           spacetime_unit_normal, phi, pi);
     910             :   auto& dt_shift =
     911             :       get<::Tags::dt<gr::Tags::Shift<DataVector, 3>>>(computation_variables);
     912             :   gh::time_deriv_of_shift(make_not_null(&dt_shift), lapse, shift,
     913             :                           inverse_spatial_metric, spacetime_unit_normal, phi,
     914             :                           pi);
     915             : 
     916             :   auto& du_null_l = get<::Tags::dt<Tags::detail::NullL>>(computation_variables);
     917             :   auto& null_l = get<Tags::detail::NullL>(computation_variables);
     918             :   null_vector_l_and_derivatives(make_not_null(&du_null_l),
     919             :                                 make_not_null(&null_l), dt_worldtube_normal,
     920             :                                 dt_lapse, dt_spacetime_metric, dt_shift, lapse,
     921             :                                 spacetime_metric, shift, worldtube_normal);
     922             : 
     923             :   // pass to the next step that is common between the 'modal' input and 'GH'
     924             :   // input strategies
     925             :   detail::create_bondi_boundary_data(
     926             :       bondi_boundary_data, make_not_null(&computation_variables),
     927             :       make_not_null(&derivative_buffers), dt_spacetime_metric, phi,
     928             :       spacetime_metric, null_l, du_null_l, cartesian_to_spherical_jacobian,
     929             :       l_max, extraction_radius);
     930             : }
     931             : 
     932             : /*!
     933             :  * \brief Process the worldtube data from modal metric components and
     934             :  * derivatives to desired Bondi quantities, placing the result in the passed
     935             :  * `Variables`.
     936             :  *
     937             :  * \details
     938             :  * The mathematics are a bit complicated for all of the coordinate
     939             :  * transformations that are necessary to obtain the Bondi gauge quantities.
     940             :  * For full mathematical details, see the documentation for functions in
     941             :  * `BoundaryData.hpp` and \cite Barkett2019uae \cite Bishop1998uk.
     942             :  *
     943             :  * This function takes as input the full set of ADM metric data and its radial
     944             :  * and time derivatives on a two-dimensional surface of constant \f$r\f$ and
     945             :  * \f$t\f$ in numerical coordinates. This data must be provided as spherical
     946             :  * harmonic coefficients in the libsharp format. This data is provided in nine
     947             :  * `Tensor`s.
     948             :  *
     949             :  * Sufficient tags to provide full worldtube boundary data at a particular
     950             :  * time are set in `bondi_boundary_data`. In particular, the set of tags in
     951             :  * `Tags::characteristic_worldtube_boundary_tags` in the provided `Variables`
     952             :  * are assigned to the worldtube boundary values associated with the input
     953             :  * metric components.
     954             :  *
     955             :  * The majority of the mathematical transformations are implemented as a set of
     956             :  * individual cascaded functions below. The details of the manipulations that
     957             :  * are performed to the input data may be found in the individual functions
     958             :  * themselves, which are called in the following order:
     959             :  * - `trigonometric_functions_on_swsh_collocation()`
     960             :  * - `cartesian_to_spherical_coordinates_and_jacobians()`
     961             :  * - `cartesian_spatial_metric_and_derivatives_from_modes()`
     962             :  * - `cartesian_shift_and_derivatives_from_modes()`
     963             :  * - `cartesian_lapse_and_derivatives_from_modes()`
     964             :  * - `gh::phi()`
     965             :  * - `gr::time_derivative_of_spacetime_metric`
     966             :  * - `gr::spacetime_metric`
     967             :  * - `generalized_harmonic_quantities()`
     968             :  * - `worldtube_normal_and_derivatives()`
     969             :  * - `null_vector_l_and_derivatives()`
     970             :  * - `null_metric_and_derivative()`
     971             :  * - `dlambda_null_metric_and_inverse()`
     972             :  * - `bondi_r()`
     973             :  * - `d_bondi_r()`
     974             :  * - `dyads()`
     975             :  * - `beta_worldtube_data()`
     976             :  * - `bondi_u_worldtube_data()`
     977             :  * - `bondi_w_worldtube_data()`
     978             :  * - `bondi_j_worldtube_data()`
     979             :  * - `dr_bondi_j()`
     980             :  * - `d2lambda_bondi_r()`
     981             :  * - `bondi_q_worldtube_data()`
     982             :  * - `bondi_h_worldtube_data()`
     983             :  * - `du_j_worldtube_data()`
     984             :  */
     985             : template <typename BoundaryTagList>
     986           1 : void create_bondi_boundary_data(
     987             :     const gsl::not_null<Variables<BoundaryTagList>*> bondi_boundary_data,
     988             :     const tnsr::ii<ComplexModalVector, 3>& spatial_metric_coefficients,
     989             :     const tnsr::ii<ComplexModalVector, 3>& dt_spatial_metric_coefficients,
     990             :     const tnsr::ii<ComplexModalVector, 3>& dr_spatial_metric_coefficients,
     991             :     const tnsr::I<ComplexModalVector, 3>& shift_coefficients,
     992             :     const tnsr::I<ComplexModalVector, 3>& dt_shift_coefficients,
     993             :     const tnsr::I<ComplexModalVector, 3>& dr_shift_coefficients,
     994             :     const Scalar<ComplexModalVector>& lapse_coefficients,
     995             :     const Scalar<ComplexModalVector>& dt_lapse_coefficients,
     996             :     const Scalar<ComplexModalVector>& dr_lapse_coefficients,
     997             :     const double extraction_radius, const size_t l_max) {
     998             :   const size_t size = Spectral::Swsh::number_of_swsh_collocation_points(l_max);
     999             : 
    1000             :   // Most allocations required for the full boundary computation are merged into
    1001             :   // a single, large Variables allocation. There remain a handful of cases in
    1002             :   // the computational functions called where an intermediate quantity that is
    1003             :   // not re-used is allocated rather than taking a buffer. These cases are
    1004             :   // marked with code comments 'Allocation'; In the future, if allocations are
    1005             :   // identified as a point to optimize, those buffers may be allocated here and
    1006             :   // passed as function arguments
    1007             :   Variables<tmpl::list<
    1008             :       Tags::detail::CosPhi, Tags::detail::CosTheta, Tags::detail::SinPhi,
    1009             :       Tags::detail::SinTheta, Tags::detail::CartesianCoordinates,
    1010             :       Tags::detail::CartesianToSphericalJacobian,
    1011             :       Tags::detail::InverseCartesianToSphericalJacobian,
    1012             :       gr::Tags::SpatialMetric<DataVector, 3>,
    1013             :       gr::Tags::InverseSpatialMetric<DataVector, 3>,
    1014             :       ::Tags::deriv<gr::Tags::SpatialMetric<DataVector, 3>, tmpl::size_t<3>,
    1015             :                     ::Frame::Inertial>,
    1016             :       ::Tags::dt<gr::Tags::SpatialMetric<DataVector, 3>>,
    1017             :       gr::Tags::Shift<DataVector, 3>,
    1018             :       ::Tags::deriv<gr::Tags::Shift<DataVector, 3>, tmpl::size_t<3>,
    1019             :                     ::Frame::Inertial>,
    1020             :       ::Tags::dt<gr::Tags::Shift<DataVector, 3>>, gr::Tags::Lapse<DataVector>,
    1021             :       ::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
    1022             :                     ::Frame::Inertial>,
    1023             :       ::Tags::dt<gr::Tags::Lapse<DataVector>>,
    1024             :       gr::Tags::SpacetimeMetric<DataVector, 3>,
    1025             :       ::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3>>,
    1026             :       gh::Tags::Phi<DataVector, 3>, Tags::detail::WorldtubeNormal,
    1027             :       ::Tags::dt<Tags::detail::WorldtubeNormal>, Tags::detail::NullL,
    1028             :       ::Tags::dt<Tags::detail::NullL>,
    1029             :       // for the detail function called at the end
    1030             :       gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>,
    1031             :       ::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>>,
    1032             :       gr::Tags::InverseSpacetimeMetric<DataVector, 3, Frame::RadialNull>,
    1033             :       Tags::detail::AngularDNullL,
    1034             :       Tags::detail::DLambda<
    1035             :           gr::Tags::SpacetimeMetric<DataVector, 3, Frame::RadialNull>>,
    1036             :       Tags::detail::DLambda<
    1037             :           gr::Tags::InverseSpacetimeMetric<DataVector, 3, Frame::RadialNull>>,
    1038             :       ::Tags::spacetime_deriv<Tags::detail::RealBondiR, tmpl::size_t<3>,
    1039             :                               Frame::RadialNull>,
    1040             :       Tags::detail::DLambda<Tags::detail::DLambda<Tags::detail::RealBondiR>>,
    1041             :       ::Tags::deriv<Tags::detail::DLambda<Tags::detail::RealBondiR>,
    1042             :                     tmpl::size_t<2>, Frame::RadialNull>>>
    1043             :       computation_variables{size};
    1044             : 
    1045             :   Variables<
    1046             :       tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
    1047             :                                       std::integral_constant<int, 0>>,
    1048             :                  ::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
    1049             :                                       std::integral_constant<int, 1>>>>
    1050             :       derivative_buffers{size};
    1051             :   auto& cos_phi = get<Tags::detail::CosPhi>(computation_variables);
    1052             :   auto& cos_theta = get<Tags::detail::CosTheta>(computation_variables);
    1053             :   auto& sin_phi = get<Tags::detail::SinPhi>(computation_variables);
    1054             :   auto& sin_theta = get<Tags::detail::SinTheta>(computation_variables);
    1055             :   trigonometric_functions_on_swsh_collocation(
    1056             :       make_not_null(&cos_phi), make_not_null(&cos_theta),
    1057             :       make_not_null(&sin_phi), make_not_null(&sin_theta), l_max);
    1058             : 
    1059             :   // NOTE: to handle the singular values of polar coordinates, the phi
    1060             :   // components of all tensors are scaled according to their sin(theta)
    1061             :   // prefactors.
    1062             :   // so, any down-index component get<2>(A) represents 1/sin(theta) A_\phi,
    1063             :   // and any up-index component get<2>(A) represents sin(theta) A^\phi.
    1064             :   // This holds for Jacobians, and so direct application of the Jacobians
    1065             :   // brings the factors through.
    1066             :   auto& cartesian_coords =
    1067             :       get<Tags::detail::CartesianCoordinates>(computation_variables);
    1068             :   auto& cartesian_to_spherical_jacobian =
    1069             :       get<Tags::detail::CartesianToSphericalJacobian>(computation_variables);
    1070             :   auto& inverse_cartesian_to_spherical_jacobian =
    1071             :       get<Tags::detail::InverseCartesianToSphericalJacobian>(
    1072             :           computation_variables);
    1073             :   cartesian_to_spherical_coordinates_and_jacobians(
    1074             :       make_not_null(&cartesian_coords),
    1075             :       make_not_null(&cartesian_to_spherical_jacobian),
    1076             :       make_not_null(&inverse_cartesian_to_spherical_jacobian), cos_phi,
    1077             :       cos_theta, sin_phi, sin_theta, extraction_radius);
    1078             : 
    1079             :   auto& cartesian_spatial_metric =
    1080             :       get<gr::Tags::SpatialMetric<DataVector, 3>>(computation_variables);
    1081             :   auto& inverse_spatial_metric =
    1082             :       get<gr::Tags::InverseSpatialMetric<DataVector, 3>>(computation_variables);
    1083             :   auto& d_cartesian_spatial_metric =
    1084             :       get<::Tags::deriv<gr::Tags::SpatialMetric<DataVector, 3>, tmpl::size_t<3>,
    1085             :                         ::Frame::Inertial>>(computation_variables);
    1086             :   auto& dt_cartesian_spatial_metric =
    1087             :       get<::Tags::dt<gr::Tags::SpatialMetric<DataVector, 3>>>(
    1088             :           computation_variables);
    1089             :   auto& interpolation_buffer =
    1090             :       get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
    1091             :                                std::integral_constant<int, 0>>>(
    1092             :           derivative_buffers);
    1093             :   Scalar<SpinWeighted<ComplexModalVector, 0>> interpolation_modal_buffer{size};
    1094             :   auto& eth_buffer =
    1095             :       get<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
    1096             :                                std::integral_constant<int, 1>>>(
    1097             :           derivative_buffers);
    1098             :   cartesian_spatial_metric_and_derivatives_from_modes(
    1099             :       make_not_null(&cartesian_spatial_metric),
    1100             :       make_not_null(&inverse_spatial_metric),
    1101             :       make_not_null(&d_cartesian_spatial_metric),
    1102             :       make_not_null(&dt_cartesian_spatial_metric),
    1103             :       make_not_null(&interpolation_modal_buffer),
    1104             :       make_not_null(&interpolation_buffer), make_not_null(&eth_buffer),
    1105             :       spatial_metric_coefficients, dr_spatial_metric_coefficients,
    1106             :       dt_spatial_metric_coefficients, inverse_cartesian_to_spherical_jacobian,
    1107             :       l_max);
    1108             : 
    1109             :   auto& cartesian_shift =
    1110             :       get<gr::Tags::Shift<DataVector, 3>>(computation_variables);
    1111             :   auto& d_cartesian_shift =
    1112             :       get<::Tags::deriv<gr::Tags::Shift<DataVector, 3>, tmpl::size_t<3>,
    1113             :                         ::Frame::Inertial>>(computation_variables);
    1114             :   auto& dt_cartesian_shift =
    1115             :       get<::Tags::dt<gr::Tags::Shift<DataVector, 3>>>(computation_variables);
    1116             : 
    1117             :   cartesian_shift_and_derivatives_from_modes(
    1118             :       make_not_null(&cartesian_shift), make_not_null(&d_cartesian_shift),
    1119             :       make_not_null(&dt_cartesian_shift),
    1120             :       make_not_null(&interpolation_modal_buffer),
    1121             :       make_not_null(&interpolation_buffer), make_not_null(&eth_buffer),
    1122             :       shift_coefficients, dr_shift_coefficients, dt_shift_coefficients,
    1123             :       inverse_cartesian_to_spherical_jacobian, l_max);
    1124             : 
    1125             :   auto& cartesian_lapse =
    1126             :       get<gr::Tags::Lapse<DataVector>>(computation_variables);
    1127             :   auto& d_cartesian_lapse =
    1128             :       get<::Tags::deriv<gr::Tags::Lapse<DataVector>, tmpl::size_t<3>,
    1129             :                         ::Frame::Inertial>>(computation_variables);
    1130             :   auto& dt_cartesian_lapse =
    1131             :       get<::Tags::dt<gr::Tags::Lapse<DataVector>>>(computation_variables);
    1132             :   cartesian_lapse_and_derivatives_from_modes(
    1133             :       make_not_null(&cartesian_lapse), make_not_null(&d_cartesian_lapse),
    1134             :       make_not_null(&dt_cartesian_lapse),
    1135             :       make_not_null(&interpolation_modal_buffer),
    1136             :       make_not_null(&interpolation_buffer), make_not_null(&eth_buffer),
    1137             :       lapse_coefficients, dr_lapse_coefficients, dt_lapse_coefficients,
    1138             :       inverse_cartesian_to_spherical_jacobian, l_max);
    1139             : 
    1140             :   auto& phi = get<gh::Tags::Phi<DataVector, 3>>(computation_variables);
    1141             :   auto& dt_spacetime_metric =
    1142             :       get<::Tags::dt<gr::Tags::SpacetimeMetric<DataVector, 3>>>(
    1143             :           computation_variables);
    1144             :   auto& spacetime_metric =
    1145             :       get<gr::Tags::SpacetimeMetric<DataVector, 3>>(computation_variables);
    1146             :   gh::phi(make_not_null(&phi), cartesian_lapse, d_cartesian_lapse,
    1147             :           cartesian_shift, d_cartesian_shift, cartesian_spatial_metric,
    1148             :           d_cartesian_spatial_metric);
    1149             :   gr::time_derivative_of_spacetime_metric(
    1150             :       make_not_null(&dt_spacetime_metric), cartesian_lapse, dt_cartesian_lapse,
    1151             :       cartesian_shift, dt_cartesian_shift, cartesian_spatial_metric,
    1152             :       dt_cartesian_spatial_metric);
    1153             :   gr::spacetime_metric(make_not_null(&spacetime_metric), cartesian_lapse,
    1154             :                        cartesian_shift, cartesian_spatial_metric);
    1155             : 
    1156             :   auto& dt_worldtube_normal =
    1157             :       get<::Tags::dt<Tags::detail::WorldtubeNormal>>(computation_variables);
    1158             :   auto& worldtube_normal =
    1159             :       get<Tags::detail::WorldtubeNormal>(computation_variables);
    1160             :   worldtube_normal_and_derivatives(
    1161             :       make_not_null(&worldtube_normal), make_not_null(&dt_worldtube_normal),
    1162             :       cos_phi, cos_theta, spacetime_metric, dt_spacetime_metric, sin_phi,
    1163             :       sin_theta, inverse_spatial_metric);
    1164             : 
    1165             :   auto& du_null_l = get<::Tags::dt<Tags::detail::NullL>>(computation_variables);
    1166             :   auto& null_l = get<Tags::detail::NullL>(computation_variables);
    1167             :   null_vector_l_and_derivatives(
    1168             :       make_not_null(&du_null_l), make_not_null(&null_l), dt_worldtube_normal,
    1169             :       dt_cartesian_lapse, dt_spacetime_metric, dt_cartesian_shift,
    1170             :       cartesian_lapse, spacetime_metric, cartesian_shift, worldtube_normal);
    1171             : 
    1172             :   // pass to the next step that is common between the 'modal' input and 'GH'
    1173             :   // input strategies
    1174             :   detail::create_bondi_boundary_data(
    1175             :       bondi_boundary_data, make_not_null(&computation_variables),
    1176             :       make_not_null(&derivative_buffers), dt_spacetime_metric, phi,
    1177             :       spacetime_metric, null_l, du_null_l, cartesian_to_spherical_jacobian,
    1178             :       l_max, extraction_radius);
    1179             : }
    1180             : }  // namespace Cce

Generated by: LCOV version 1.14