SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/Cce - Equations.hpp Hit Total Coverage
Commit: 3c072f0ce967e2e56649d3fa12aa2a0e4fe2a42e Lines: 14 109 12.8 %
Date: 2024-04-23 20:50:18
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 <type_traits>
       7             : 
       8             : #include "DataStructures/SpinWeighted.hpp"  // IWYU pragma: keep
       9             : #include "DataStructures/Tags.hpp"
      10             : #include "DataStructures/Tags/TempTensor.hpp"
      11             : #include "DataStructures/Tensor/Tensor.hpp"
      12             : #include "Evolution/Systems/Cce/Tags.hpp"  // IWYU pragma: keep
      13             : #include "NumericalAlgorithms/Spectral/SwshTags.hpp"
      14             : #include "Utilities/Gsl.hpp"
      15             : #include "Utilities/TMPL.hpp"
      16             : 
      17             : // IWYU pragma: no_forward_declare Cce::Tags::BondiBeta
      18             : // IWYU pragma: no_forward_declare Cce::Tags::DuRDividedByR
      19             : // IWYU pragma: no_forward_declare Cce::Tags::EthRDividedByR
      20             : // IWYU pragma: no_forward_declare Cce::Tags::Exp2Beta
      21             : // IWYU pragma: no_forward_declare Cce::Tags::BondiH
      22             : // IWYU pragma: no_forward_declare Cce::Tags::BondiJ
      23             : // IWYU pragma: no_forward_declare Cce::Tags::BondiJbar
      24             : // IWYU pragma: no_forward_declare Cce::Tags::JbarQMinus2EthBeta
      25             : // IWYU pragma: no_forward_declare Cce::Tags::BondiK
      26             : // IWYU pragma: no_forward_declare Cce::Tags::OneMinusY
      27             : // IWYU pragma: no_forward_declare Cce::Tags::BondiQ
      28             : // IWYU pragma: no_forward_declare Cce::Tags::BondiR
      29             : // IWYU pragma: no_forward_declare Cce::Tags::BondiU
      30             : // IWYU pragma: no_forward_declare Cce::Tags::BondiUbar
      31             : // IWYU pragma: no_forward_declare Cce::Tags::BondiW
      32             : // IWYU pragma: no_forward_declare ::Tags::Multiplies
      33             : // IWYU pragma: no_forward_declare Cce::Tags::Dy
      34             : // IWYU pragma: no_forward_declare Cce::Tags::Integrand
      35             : // IWYU pragma: no_forward_declare Cce::Tags::LinearFactor
      36             : // IWYU pragma: no_forward_declare Cce::Tags::LinearFactorForConjugate
      37             : // IWYU pragma: no_forward_declare Cce::Tags::PoleOfIntegrand
      38             : // IWYU pragma: no_forward_declare Cce::Tags::RegularIntegrand
      39             : // IWYU pragma: no_forward_declare Spectral::Swsh::Tags::Eth
      40             : // IWYU pragma: no_forward_declare Spectral::Swsh::Tags::EthEth
      41             : // IWYU pragma: no_forward_declare Spectral::Swsh::Tags::EthEthbar
      42             : // IWYU pragma: no_forward_declare Spectral::Swsh::Tags::Ethbar
      43             : // IWYU pragma: no_forward_declare Spectral::Swsh::Tags::EthbarEthbar
      44             : // IWYU pragma: no_forward_declare Spectral::Swsh::Tags::Derivative
      45             : // IWYU pragma: no_forward_declare Tags::TempTensor
      46             : // IWYU pragma: no_forward_declare Tags::SpinWeighted
      47             : // IWYU pragma: no_forward_declare SpinWeighted
      48             : // IWYU pragma: no_forward_declare Tensor
      49             : 
      50             : /// \cond
      51             : class ComplexDataVector;
      52             : /// \endcond
      53             : 
      54             : namespace Cce {
      55             : 
      56             : namespace detail {
      57             : template <typename BondiVariable>
      58             : struct integrand_terms_to_compute_for_bondi_variable_impl;
      59             : 
      60             : // template specializations for the individual tags
      61             : template <>
      62             : struct integrand_terms_to_compute_for_bondi_variable_impl<Tags::BondiBeta> {
      63             :   using type = tmpl::list<Tags::Integrand<Tags::BondiBeta>>;
      64             : };
      65             : template <>
      66             : struct integrand_terms_to_compute_for_bondi_variable_impl<Tags::BondiQ> {
      67             :   using type = tmpl::list<Tags::PoleOfIntegrand<Tags::BondiQ>,
      68             :                           Tags::RegularIntegrand<Tags::BondiQ>>;
      69             : };
      70             : template <>
      71             : struct integrand_terms_to_compute_for_bondi_variable_impl<Tags::BondiU> {
      72             :   using type = tmpl::list<Tags::Integrand<Tags::BondiU>>;
      73             : };
      74             : template <>
      75             : struct integrand_terms_to_compute_for_bondi_variable_impl<Tags::BondiW> {
      76             :   using type = tmpl::list<Tags::PoleOfIntegrand<Tags::BondiW>,
      77             :                           Tags::RegularIntegrand<Tags::BondiW>>;
      78             : };
      79             : template <>
      80             : struct integrand_terms_to_compute_for_bondi_variable_impl<Tags::BondiH> {
      81             :   using type = tmpl::list<Tags::PoleOfIntegrand<Tags::BondiH>,
      82             :                           Tags::RegularIntegrand<Tags::BondiH>,
      83             :                           Tags::LinearFactor<Tags::BondiH>,
      84             :                           Tags::LinearFactorForConjugate<Tags::BondiH>>;
      85             : };
      86             : 
      87             : template <>
      88             : struct integrand_terms_to_compute_for_bondi_variable_impl<Tags::KleinGordonPi> {
      89             :   using type = tmpl::list<Tags::PoleOfIntegrand<Tags::KleinGordonPi>,
      90             :                           Tags::RegularIntegrand<Tags::KleinGordonPi>>;
      91             : };
      92             : 
      93             : void klein_gordon_rhs_npsi1(
      94             :     gsl::not_null<SpinWeighted<ComplexDataVector, 0>*> result,
      95             :     const SpinWeighted<ComplexDataVector, 1>& eth_beta,
      96             :     const SpinWeighted<ComplexDataVector, 1>& eth_kg_psi,
      97             :     const SpinWeighted<ComplexDataVector, 0>& eth_ethbar_kg_psi,
      98             :     const SpinWeighted<ComplexDataVector, 0>& bondi_k);
      99             : 
     100             : void klein_gordon_rhs_npsi2(
     101             :     gsl::not_null<SpinWeighted<ComplexDataVector, 0>*> result,
     102             :     const SpinWeighted<ComplexDataVector, 2>& j,
     103             :     const SpinWeighted<ComplexDataVector, 2>& eth_eth_kg_psi,
     104             :     const SpinWeighted<ComplexDataVector, 1>& eth_beta,
     105             :     const SpinWeighted<ComplexDataVector, 1>& eth_kg_psi,
     106             :     const SpinWeighted<ComplexDataVector, 1>& ethbar_j);
     107             : 
     108             : void klein_gordon_rhs_npsi3(
     109             :     gsl::not_null<SpinWeighted<ComplexDataVector, 0>*> result,
     110             :     const SpinWeighted<ComplexDataVector, 1>& eth_k,
     111             :     const SpinWeighted<ComplexDataVector, 1>& eth_kg_psi);
     112             : 
     113             : void klein_gordon_rhs_npsi4_divided_by_one_minues_y_squared(
     114             :     gsl::not_null<SpinWeighted<ComplexDataVector, 0>*> result,
     115             :     const SpinWeighted<ComplexDataVector, 1>& eth_kg_psi,
     116             :     const SpinWeighted<ComplexDataVector, 1>& dy_bondi_u,
     117             :     const SpinWeighted<ComplexDataVector, 0>& ethbar_u,
     118             :     const SpinWeighted<ComplexDataVector, 1>& bondi_u,
     119             :     const SpinWeighted<ComplexDataVector, 1>& eth_dy_kg_psi,
     120             :     const SpinWeighted<ComplexDataVector, 0>& dy_kg_psi,
     121             :     const SpinWeighted<ComplexDataVector, 0>& bondi_r,
     122             :     const SpinWeighted<ComplexDataVector, 1>& eth_r_divided_by_r);
     123             : 
     124             : void klein_gordon_rhs_tau(
     125             :     gsl::not_null<SpinWeighted<ComplexDataVector, 0>*> result,
     126             :     const SpinWeighted<ComplexDataVector, 0>& one_minus_y,
     127             :     const SpinWeighted<ComplexDataVector, 0>& dy_w,
     128             :     const SpinWeighted<ComplexDataVector, 0>& dy_kg_psi,
     129             :     const SpinWeighted<ComplexDataVector, 0>& dy_dy_kg_psi,
     130             :     const SpinWeighted<ComplexDataVector, 0>& bondi_w,
     131             :     const SpinWeighted<ComplexDataVector, 0>& bondi_r);
     132             : }  // namespace detail
     133             : 
     134             : /// \brief A struct for providing a `tmpl::list` of integrand tags that need to
     135             : /// be computed before integration can proceed for a given Bondi variable tag.
     136             : template <typename BondiVariable>
     137           1 : using integrand_terms_to_compute_for_bondi_variable =
     138             :     typename detail::integrand_terms_to_compute_for_bondi_variable_impl<
     139             :         BondiVariable>::type;
     140             : 
     141             : /*!
     142             :  * \brief Computes one of the inputs for the integration of one of the
     143             :  * Characteristic hypersurface equations.
     144             :  *
     145             :  * \details The template argument must be one of the integrand-related prefix
     146             :  * tags templated on a Bondi quantity tag for which that integrand is required.
     147             :  * The relevant prefix tags are `Tags::Integrand`, `Tags::PoleOfIntegrand`,
     148             :  * `Tags::RegularIntegrand`, `Tags::LinearFactor`, and
     149             :  * `Tags::LinearFactorForConjugate`. The Bondi quantity tags that these tags may
     150             :  * wrap are `Tags::BondiBeta`, `Tags::BondiQ`, `Tags::BondiU`, `Tags::BondiW`,
     151             :  * and `Tags::BondiH`.
     152             :  *
     153             :  * The integrand terms which may be computed for a given Bondi variable are
     154             :  * enumerated in the type alias `integrand_terms_to_compute_for_bondi_variable`,
     155             :  * which takes as a single template argument the tag for which integrand terms
     156             :  * would be computed, and is a `tmpl::list` of the integrand terms needed.
     157             :  *
     158             :  * The resulting quantity is returned by `not_null` pointer, and the required
     159             :  * argument tags are given in `return_tags` and `argument_tags` type aliases,
     160             :  * where the `return_tags` are passed by `not_null` pointer (so include
     161             :  * temporary buffers) and the `argument_tags` are passed by const reference.
     162             :  *
     163             :  * Additional mathematical details for each of the computations can be found in
     164             :  * the template specializations of this struct. All of the specializations have
     165             :  * a static `apply` function which takes as arguments
     166             :  * `SpinWeighted<ComplexDataVector, N>`s for each of the Bondi arguments.
     167             :  */
     168             : template <typename IntegrandTag>
     169           1 : struct ComputeBondiIntegrand;
     170             : 
     171             : /*!
     172             :  * \brief Computes the integrand (right-hand side) of the equation which
     173             :  * determines the radial (y) dependence of the Bondi quantity \f$\beta\f$.
     174             :  *
     175             :  * \details The quantity \f$\beta\f$ is defined via the Bondi form of the
     176             :  * metric:
     177             :  * \f[
     178             :  * ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2
     179             :  * e^{2 \beta} du dr - 2 r^2 h_{AB} U^B du dx^A + r^2 h_{A B} dx^A dx^B. \f]
     180             :  * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical
     181             :  * angular dyad \f$q^A\f$:
     182             :  * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A
     183             :  * \bar{q}^B.\f]
     184             :  * See \cite Bishop1997ik \cite Handmer2014qha for full details.
     185             :  *
     186             :  * We write the equations of motion in the compactified coordinate \f$ y \equiv
     187             :  * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the
     188             :  * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of
     189             :  * the worldtube. The equation which determines \f$\beta\f$ on a surface of
     190             :  * constant \f$u\f$ given \f$J\f$ on the same surface is
     191             :  * \f[\partial_y (\beta) =
     192             :  * \frac{1}{8} (-1 + y) \left(\partial_y (J) \partial_y(\bar{J})
     193             :  * - \frac{(\partial_y (J \bar{J}))^2}{4 K^2}\right). \f]
     194             :  */
     195             : template <>
     196           1 : struct ComputeBondiIntegrand<Tags::Integrand<Tags::BondiBeta>> {
     197             :  public:
     198           0 :   using pre_swsh_derivative_tags =
     199             :       tmpl::list<Tags::Dy<Tags::BondiJ>, Tags::BondiJ>;
     200           0 :   using swsh_derivative_tags = tmpl::list<>;
     201           0 :   using integration_independent_tags = tmpl::list<Tags::OneMinusY>;
     202           0 :   using temporary_tags = tmpl::list<>;
     203             : 
     204           0 :   using return_tags = tmpl::append<tmpl::list<Tags::Integrand<Tags::BondiBeta>>,
     205             :                                    temporary_tags>;
     206           0 :   using argument_tags =
     207             :       tmpl::append<pre_swsh_derivative_tags, swsh_derivative_tags,
     208             :                    integration_independent_tags>;
     209             : 
     210             :   template <typename... Args>
     211           0 :   static void apply(
     212             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
     213             :           integrand_for_beta,
     214             :       const Args&... args) {
     215             :     apply_impl(make_not_null(&get(*integrand_for_beta)), get(args)...);
     216             :   }
     217             : 
     218             :  private:
     219           0 :   static void apply_impl(
     220             :       gsl::not_null<SpinWeighted<ComplexDataVector, 0>*> integrand_for_beta,
     221             :       const SpinWeighted<ComplexDataVector, 2>& dy_j,
     222             :       const SpinWeighted<ComplexDataVector, 2>& j,
     223             :       const SpinWeighted<ComplexDataVector, 0>& one_minus_y);
     224             : };
     225             : 
     226             : /*!
     227             :  * \brief Computes the pole part of the integrand (right-hand side) of the
     228             :  * equation which determines the radial (y) dependence of the Bondi quantity
     229             :  * \f$Q\f$.
     230             :  *
     231             :  * \details The quantity \f$Q\f$ is defined via the Bondi form of the metric:
     232             :  * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2
     233             :  * e^{2 \beta} du dr - 2 r^2 h_{AB} U^B du dx^A + r^2 h_{A B} dx^A dx^B. \f]
     234             :  * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical
     235             :  * angular dyad \f$q^A\f$:
     236             :  * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B,\f]
     237             :  * and \f$Q\f$ is defined as a supplemental variable for radial integration of
     238             :  * \f$U\f$:
     239             :  * \f[ Q_A = r^2 e^{-2\beta} h_{AB} \partial_r U^B\f]
     240             :  * and \f$Q = Q_A q^A\f$. See \cite Bishop1997ik \cite Handmer2014qha for full
     241             :  * details.
     242             :  *
     243             :  * We write the equations of motion in the compactified coordinate \f$ y \equiv
     244             :  * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the
     245             :  * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of
     246             :  * the worldtube. The equation which determines \f$Q\f$ on a surface of constant
     247             :  * \f$u\f$ given \f$J\f$ and \f$\beta\f$ on the same surface is written as
     248             :  * \f[(1 - y) \partial_y Q + 2 Q = A_Q + (1 - y) B_Q.\f]
     249             :  * We refer to \f$A_Q\f$ as the "pole part" of the integrand and \f$B_Q\f$
     250             :  * as the "regular part". The pole part is computed by this function, and has
     251             :  * the expression
     252             :  * \f[A_Q = -4 \eth \beta.\f]
     253             :  */
     254             : template <>
     255           1 : struct ComputeBondiIntegrand<Tags::PoleOfIntegrand<Tags::BondiQ>> {
     256             :  public:
     257           0 :   using pre_swsh_derivative_tags = tmpl::list<>;
     258           0 :   using swsh_derivative_tags =
     259             :       tmpl::list<Spectral::Swsh::Tags::Derivative<Tags::BondiBeta,
     260             :                                                   Spectral::Swsh::Tags::Eth>>;
     261           0 :   using integration_independent_tags = tmpl::list<>;
     262           0 :   using temporary_tags = tmpl::list<>;
     263             : 
     264           0 :   using return_tags =
     265             :       tmpl::append<tmpl::list<Tags::PoleOfIntegrand<Tags::BondiQ>>,
     266             :                    temporary_tags>;
     267           0 :   using argument_tags =
     268             :       tmpl::append<pre_swsh_derivative_tags, swsh_derivative_tags,
     269             :                    integration_independent_tags>;
     270             : 
     271             :   template <typename... Args>
     272           0 :   static void apply(
     273             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*>
     274             :           pole_of_integrand_for_q,
     275             :       const Args&... args) {
     276             :     apply_impl(make_not_null(&get(*pole_of_integrand_for_q)), get(args)...);
     277             :   }
     278             : 
     279             :  private:
     280           0 :   static void apply_impl(gsl::not_null<SpinWeighted<ComplexDataVector, 1>*>
     281             :                              pole_of_integrand_for_q,
     282             :                          const SpinWeighted<ComplexDataVector, 1>& eth_beta);
     283             : };
     284             : 
     285             : /*!
     286             :  * \brief Computes the regular part of the integrand (right-hand side) of the
     287             :  * equation which determines the radial (y) dependence of the Bondi quantity
     288             :  * \f$Q\f$.
     289             :  *
     290             :  * \details The quantity \f$Q\f$ is defined via the Bondi form of the metric:
     291             :  * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2
     292             :  * e^{2 \beta} du dr - 2 r^2 h_{AB} U^B du dx^A + r^2 h_{A B} dx^A dx^B. \f]
     293             :  * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical
     294             :  * angular dyad \f$q^A\f$:
     295             :  * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B,\f]
     296             :  * and \f$Q\f$ is defined as a supplemental variable for radial integration of
     297             :  * \f$U\f$:
     298             :  * \f[ Q_A = r^2 e^{-2\beta} h_{AB} \partial_r U^B\f]
     299             :  * and \f$Q = Q_A q^A\f$. See \cite Bishop1997ik \cite Handmer2014qha for
     300             :  * full details.
     301             :  *
     302             :  * We write the equations of motion in the compactified coordinate \f$ y \equiv
     303             :  * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the
     304             :  * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of
     305             :  * the worldtube. The equation which determines \f$Q\f$ on a surface of constant
     306             :  * \f$u\f$ given \f$J\f$ and \f$\beta\f$ on the same surface is written as
     307             :  * \f[(1 - y) \partial_y Q + 2 Q = A_Q + (1 - y) B_Q. \f]
     308             :  * We refer to \f$A_Q\f$ as the "pole part" of the integrand and \f$B_Q\f$ as
     309             :  * the "regular part". The regular part is computed by this function, and has
     310             :  * the expression
     311             :  * \f[  B_Q = - \left(2 \mathcal{A}_Q + \frac{2
     312             :  * \bar{\mathcal{A}_Q} J}{K} -  2 \partial_y (\eth (\beta)) +
     313             :  * \frac{\partial_y (\bar{\eth} (J))}{K}\right), \f]
     314             :  * where
     315             :  * \f[ \mathcal{A}_Q = - \tfrac{1}{4} \eth (\bar{J} \partial_y (J)) +
     316             :  * \tfrac{1}{4} J \partial_y (\eth (\bar{J})) -  \tfrac{1}{4} \eth (\bar{J})
     317             :  * \partial_y (J) + \frac{\eth (J \bar{J}) \partial_y (J \bar{J})}{8 K^2} -
     318             :  * \frac{\bar{J} \eth (R) \partial_y (J)}{4 R}. \f].
     319             :  */
     320             : template <>
     321           1 : struct ComputeBondiIntegrand<Tags::RegularIntegrand<Tags::BondiQ>> {
     322             :  public:
     323           0 :   using pre_swsh_derivative_tags =
     324             :       tmpl::list<Tags::Dy<Tags::BondiBeta>, Tags::Dy<Tags::BondiJ>,
     325             :                  Tags::BondiJ>;
     326           0 :   using swsh_derivative_tags = tmpl::list<
     327             :       Spectral::Swsh::Tags::Derivative<Tags::Dy<Tags::BondiBeta>,
     328             :                                        Spectral::Swsh::Tags::Eth>,
     329             :       Spectral::Swsh::Tags::Derivative<
     330             :           ::Tags::Multiplies<Tags::BondiJ, Tags::BondiJbar>,
     331             :           Spectral::Swsh::Tags::Eth>,
     332             :       Spectral::Swsh::Tags::Derivative<
     333             :           ::Tags::Multiplies<Tags::BondiJbar, Tags::Dy<Tags::BondiJ>>,
     334             :           Spectral::Swsh::Tags::Eth>,
     335             :       Spectral::Swsh::Tags::Derivative<Tags::Dy<Tags::BondiJ>,
     336             :                                        Spectral::Swsh::Tags::Ethbar>,
     337             :       Spectral::Swsh::Tags::Derivative<Tags::BondiJ,
     338             :                                        Spectral::Swsh::Tags::Ethbar>>;
     339           0 :   using integration_independent_tags =
     340             :       tmpl::list<Tags::EthRDividedByR, Tags::BondiK>;
     341           0 :   using temporary_tags =
     342             :       tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
     343             :                                       std::integral_constant<int, 1>>>;
     344             : 
     345           0 :   using return_tags =
     346             :       tmpl::append<tmpl::list<Tags::RegularIntegrand<Tags::BondiQ>>,
     347             :                    temporary_tags>;
     348           0 :   using argument_tags =
     349             :       tmpl::append<pre_swsh_derivative_tags, swsh_derivative_tags,
     350             :                    integration_independent_tags>;
     351             : 
     352             :   template <typename... Args>
     353           0 :   static void apply(
     354             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*>
     355             :           regular_integrand_for_q,
     356             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*>
     357             :           script_aq,
     358             :       const Args&... args) {
     359             :     apply_impl(make_not_null(&get(*regular_integrand_for_q)),
     360             :                make_not_null(&get(*script_aq)), get(args)...);
     361             :   }
     362             : 
     363             :  private:
     364           0 :   static void apply_impl(
     365             :       gsl::not_null<SpinWeighted<ComplexDataVector, 1>*>
     366             :           regular_integrand_for_q,
     367             :       gsl::not_null<SpinWeighted<ComplexDataVector, 1>*> script_aq,
     368             :       const SpinWeighted<ComplexDataVector, 0>& dy_beta,
     369             :       const SpinWeighted<ComplexDataVector, 2>& dy_j,
     370             :       const SpinWeighted<ComplexDataVector, 2>& j,
     371             :       const SpinWeighted<ComplexDataVector, 1>& eth_dy_beta,
     372             :       const SpinWeighted<ComplexDataVector, 1>& eth_j_jbar,
     373             :       const SpinWeighted<ComplexDataVector, 1>& eth_jbar_dy_j,
     374             :       const SpinWeighted<ComplexDataVector, 1>& ethbar_dy_j,
     375             :       const SpinWeighted<ComplexDataVector, 1>& ethbar_j,
     376             :       const SpinWeighted<ComplexDataVector, 1>& eth_r_divided_by_r,
     377             :       const SpinWeighted<ComplexDataVector, 0>& k);
     378             : };
     379             : 
     380             : /*!
     381             :  * \brief Computes the integrand (right-hand side) of the equation which
     382             :  * determines the radial (y) dependence of the Bondi quantity \f$U\f$.
     383             :  *
     384             :  * \details The quantity \f$U\f$ is defined via the Bondi form of the metric:
     385             :  * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2
     386             :  * e^{2 \beta} du dr - 2 r^2 h_{AB} U^B du dx^A + r^2 h_{A B} dx^A dx^B. \f]
     387             :  * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical
     388             :  * angular dyad \f$q^A\f$:
     389             :  * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B,\f]
     390             :  * and \f$Q\f$ is defined as a supplemental variable for radial integration of
     391             :  * \f$U\f$:
     392             :  * \f[ Q_A = r^2 e^{-2\beta} h_{AB} \partial_r U^B\f]
     393             :  * and \f$U = U_A q^A\f$. See \cite Bishop1997ik \cite Handmer2014qha for full
     394             :  * details.
     395             :  *
     396             :  * We write the equations of motion in the compactified coordinate \f$ y \equiv
     397             :  * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the
     398             :  * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of
     399             :  * the worldtube. The equation which determines \f$U\f$ on a surface of constant
     400             :  * \f$u\f$ given \f$J\f$, \f$\beta\f$, and \f$Q\f$ on the same surface is
     401             :  * written as
     402             :  * \f[\partial_y U = \frac{e^{2\beta}}{2 R} (K Q - J \bar{Q}). \f]
     403             :  */
     404             : template <>
     405           1 : struct ComputeBondiIntegrand<Tags::Integrand<Tags::BondiU>> {
     406             :  public:
     407           0 :   using pre_swsh_derivative_tags =
     408             :       tmpl::list<Tags::Exp2Beta, Tags::BondiJ, Tags::BondiQ>;
     409           0 :   using swsh_derivative_tags = tmpl::list<>;
     410           0 :   using integration_independent_tags = tmpl::list<Tags::BondiK, Tags::BondiR>;
     411           0 :   using temporary_tags = tmpl::list<>;
     412             : 
     413           0 :   using return_tags =
     414             :       tmpl::append<tmpl::list<Tags::Integrand<Tags::BondiU>>, temporary_tags>;
     415           0 :   using argument_tags =
     416             :       tmpl::append<pre_swsh_derivative_tags, swsh_derivative_tags,
     417             :                    integration_independent_tags>;
     418             : 
     419             :   template <typename... Args>
     420           0 :   static void apply(
     421             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*>
     422             :           regular_integrand_for_u,
     423             :       const Args&... args) {
     424             :     apply_impl(make_not_null(&get(*regular_integrand_for_u)), get(args)...);
     425             :   }
     426             : 
     427             :  private:
     428           0 :   static void apply_impl(gsl::not_null<SpinWeighted<ComplexDataVector, 1>*>
     429             :                              regular_integrand_for_u,
     430             :                          const SpinWeighted<ComplexDataVector, 0>& exp_2_beta,
     431             :                          const SpinWeighted<ComplexDataVector, 2>& j,
     432             :                          const SpinWeighted<ComplexDataVector, 1>& q,
     433             :                          const SpinWeighted<ComplexDataVector, 0>& k,
     434             :                          const SpinWeighted<ComplexDataVector, 0>& r);
     435             : };
     436             : 
     437             : /*!
     438             :  * \brief Computes the pole part of the integrand (right-hand side) of the
     439             :  * equation which determines the radial (y) dependence of the Bondi quantity
     440             :  * \f$W\f$.
     441             :  *
     442             :  * \details The quantity \f$W\f$ is defined via the Bondi form of the metric:
     443             :  * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2
     444             :  * e^{2 \beta} du dr - 2 r^2 h_{AB} U^B du dx^A + r^2 h_{A B} dx^A dx^B. \f]
     445             :  * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical
     446             :  * angular dyad \f$q^A\f$:
     447             :  * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B.\f]
     448             :  * See \cite Bishop1997ik \cite Handmer2014qha for full details.
     449             :  *
     450             :  * We write the equations of motion in the compactified coordinate \f$ y \equiv
     451             :  * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the
     452             :  * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of
     453             :  * the worldtube. The equation which determines \f$W\f$ on a surface of constant
     454             :  * \f$u\f$ given \f$J\f$,\f$\beta\f$, \f$Q\f$, and \f$U\f$ on the same surface
     455             :  * is written as
     456             :  * \f[(1 - y) \partial_y W + 2 W = A_W + (1 - y) B_W.\f] We refer
     457             :  * to \f$A_W\f$ as the "pole part" of the integrand and \f$B_W\f$ as the
     458             :  * "regular part". The pole part is computed by this function, and has the
     459             :  * expression
     460             :  * \f[A_W = \eth (\bar{U}) + \bar{\eth} (U).\f]
     461             :  */
     462             : template <>
     463           1 : struct ComputeBondiIntegrand<Tags::PoleOfIntegrand<Tags::BondiW>> {
     464             :  public:
     465           0 :   using pre_swsh_derivative_tags = tmpl::list<>;
     466           0 :   using swsh_derivative_tags = tmpl::list<Spectral::Swsh::Tags::Derivative<
     467             :       Tags::BondiU, Spectral::Swsh::Tags::Ethbar>>;
     468           0 :   using integration_independent_tags = tmpl::list<>;
     469           0 :   using temporary_tags = tmpl::list<>;
     470             : 
     471           0 :   using return_tags =
     472             :       tmpl::append<tmpl::list<Tags::PoleOfIntegrand<Tags::BondiW>>,
     473             :                    temporary_tags>;
     474           0 :   using argument_tags =
     475             :       tmpl::append<pre_swsh_derivative_tags, swsh_derivative_tags,
     476             :                    integration_independent_tags>;
     477             : 
     478             :   template <typename... Args>
     479           0 :   static void apply(
     480             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
     481             :           pole_of_integrand_for_w,
     482             :       const Args&... args) {
     483             :     apply_impl(make_not_null(&get(*pole_of_integrand_for_w)), get(args)...);
     484             :   }
     485             : 
     486             :  private:
     487           0 :   static void apply_impl(gsl::not_null<SpinWeighted<ComplexDataVector, 0>*>
     488             :                              pole_of_integrand_for_w,
     489             :                          const SpinWeighted<ComplexDataVector, 0>& ethbar_u);
     490             : };
     491             : 
     492             : /*!
     493             :  * \brief Computes the regular part of the integrand (right-hand side) of the
     494             :  * equation which determines the radial (y) dependence of the Bondi quantity
     495             :  * \f$W\f$.
     496             :  *
     497             :  * \details The quantity \f$W\f$ is defined via the Bondi form of the metric:
     498             :  * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2
     499             :  * e^{2 \beta} du dr - 2 r^2 h_{AB} U^B du dx^A + r^2 h_{A B} dx^A dx^B. \f]
     500             :  * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical
     501             :  * angular dyad \f$q^A\f$:
     502             :  * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B,\f]
     503             :  * See \cite Bishop1997ik \cite Handmer2014qha for full details.
     504             :  *
     505             :  * We write the equations of motion in the compactified coordinate \f$ y \equiv
     506             :  * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the
     507             :  * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of
     508             :  * the worldtube. The equation which determines \f$W\f$ on a surface of constant
     509             :  * \f$u\f$ given \f$J\f$, \f$\beta\f$, \f$Q\f$, \f$U\f$ on the same surface is
     510             :  * written as
     511             :  * \f[(1 - y) \partial_y W + 2 W = A_W + (1 - y) B_W. \f]
     512             :  * We refer to \f$A_W\f$ as the "pole part" of the integrand and \f$B_W\f$ as
     513             :  * the "regular part". The regular part is computed by this function, and has
     514             :  * the expression
     515             :  * \f[  B_W = \tfrac{1}{4} \partial_y (\eth (\bar{U})) + \tfrac{1}{4} \partial_y
     516             :  * (\bar{\eth} (U)) -  \frac{1}{2 R} + \frac{e^{2 \beta} (\mathcal{A}_W +
     517             :  * \bar{\mathcal{A}_W})}{4 R}, \f]
     518             :  * where
     519             :  * \f{align*}
     520             :  * \mathcal{A}_W =& - \eth (\beta) \eth (\bar{J}) + \tfrac{1}{2} \bar{\eth}
     521             :  * (\bar{\eth} (J)) + 2 \bar{\eth} (\beta) \bar{\eth} (J) + (\bar{\eth}
     522             :  * (\beta))^2 J + \bar{\eth}
     523             :  * (\bar{\eth} (\beta)) J + \frac{\eth (J \bar{J}) \bar{\eth} (J \bar{J})}{8
     524             :  * K^3} + \frac{1}{2 K} -  \frac{\eth (\bar{\eth} (J \bar{J}))}{8 K} -
     525             :  * \frac{\eth (J
     526             :  * \bar{J}) \bar{\eth} (\beta)}{2 K} \nonumber \\
     527             :  * &-  \frac{\eth (\bar{J}) \bar{\eth} (J)}{4 K} -  \frac{\eth (\bar{\eth} (J))
     528             :  * \bar{J}}{4 K} + \tfrac{1}{2} K  -  \eth (\bar{\eth} (\beta)) K -  \eth
     529             :  * (\beta) \bar{\eth} (\beta) K + \tfrac{1}{4} (- K Q \bar{Q} + J \bar{Q}^2).
     530             :  * \f}
     531             :  */
     532             : template <>
     533           1 : struct ComputeBondiIntegrand<Tags::RegularIntegrand<Tags::BondiW>> {
     534             :  public:
     535           0 :   using pre_swsh_derivative_tags =
     536             :       tmpl::list<Tags::Dy<Tags::BondiU>, Tags::Exp2Beta, Tags::BondiJ,
     537             :                  Tags::BondiQ>;
     538           0 :   using swsh_derivative_tags = tmpl::list<
     539             :       Spectral::Swsh::Tags::Derivative<Tags::BondiBeta,
     540             :                                        Spectral::Swsh::Tags::Eth>,
     541             :       Spectral::Swsh::Tags::Derivative<Tags::BondiBeta,
     542             :                                        Spectral::Swsh::Tags::EthEth>,
     543             :       Spectral::Swsh::Tags::Derivative<Tags::BondiBeta,
     544             :                                        Spectral::Swsh::Tags::EthEthbar>,
     545             :       Spectral::Swsh::Tags::Derivative<
     546             :           Spectral::Swsh::Tags::Derivative<Tags::BondiJ,
     547             :                                            Spectral::Swsh::Tags::Ethbar>,
     548             :           Spectral::Swsh::Tags::Eth>,
     549             :       Spectral::Swsh::Tags::Derivative<
     550             :           ::Tags::Multiplies<Tags::BondiJ, Tags::BondiJbar>,
     551             :           Spectral::Swsh::Tags::EthEthbar>,
     552             :       Spectral::Swsh::Tags::Derivative<
     553             :           ::Tags::Multiplies<Tags::BondiJ, Tags::BondiJbar>,
     554             :           Spectral::Swsh::Tags::Eth>,
     555             :       Spectral::Swsh::Tags::Derivative<Tags::Dy<Tags::BondiU>,
     556             :                                        Spectral::Swsh::Tags::Ethbar>,
     557             :       Spectral::Swsh::Tags::Derivative<Tags::BondiJ,
     558             :                                        Spectral::Swsh::Tags::EthbarEthbar>,
     559             :       Spectral::Swsh::Tags::Derivative<Tags::BondiJ,
     560             :                                        Spectral::Swsh::Tags::Ethbar>>;
     561           0 :   using integration_independent_tags =
     562             :       tmpl::list<Tags::EthRDividedByR, Tags::BondiK, Tags::BondiR>;
     563           0 :   using temporary_tags =
     564             :       tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
     565             :                                       std::integral_constant<int, 0>>>;
     566             : 
     567           0 :   using return_tags =
     568             :       tmpl::append<tmpl::list<Tags::RegularIntegrand<Tags::BondiW>>,
     569             :                    temporary_tags>;
     570           0 :   using argument_tags =
     571             :       tmpl::append<pre_swsh_derivative_tags, swsh_derivative_tags,
     572             :                    integration_independent_tags>;
     573             : 
     574             :   template <typename... Args>
     575           0 :   static void apply(
     576             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
     577             :           regular_integrand_for_w,
     578             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
     579             :           script_av,
     580             :       const Args&... args) {
     581             :     apply_impl(make_not_null(&get(*regular_integrand_for_w)),
     582             :                make_not_null(&get(*script_av)), get(args)...);
     583             :   }
     584             : 
     585             :  private:
     586           0 :   static void apply_impl(
     587             :       gsl::not_null<SpinWeighted<ComplexDataVector, 0>*>
     588             :           regular_integrand_for_w,
     589             :       gsl::not_null<SpinWeighted<ComplexDataVector, 0>*> script_av,
     590             :       const SpinWeighted<ComplexDataVector, 1>& dy_u,
     591             :       const SpinWeighted<ComplexDataVector, 0>& exp_2_beta,
     592             :       const SpinWeighted<ComplexDataVector, 2>& j,
     593             :       const SpinWeighted<ComplexDataVector, 1>& q,
     594             :       const SpinWeighted<ComplexDataVector, 1>& eth_beta,
     595             :       const SpinWeighted<ComplexDataVector, 2>& eth_eth_beta,
     596             :       const SpinWeighted<ComplexDataVector, 0>& eth_ethbar_beta,
     597             :       const SpinWeighted<ComplexDataVector, 2>& eth_ethbar_j,
     598             :       const SpinWeighted<ComplexDataVector, 0>& eth_ethbar_j_jbar,
     599             :       const SpinWeighted<ComplexDataVector, 1>& eth_j_jbar,
     600             :       const SpinWeighted<ComplexDataVector, 0>& ethbar_dy_u,
     601             :       const SpinWeighted<ComplexDataVector, 0>& ethbar_ethbar_j,
     602             :       const SpinWeighted<ComplexDataVector, 1>& ethbar_j,
     603             :       const SpinWeighted<ComplexDataVector, 1>& eth_r_divided_by_r,
     604             :       const SpinWeighted<ComplexDataVector, 0>& k,
     605             :       const SpinWeighted<ComplexDataVector, 0>& r);
     606             : };
     607             : 
     608             : /*!
     609             :  * \brief Computes the pole part of the integrand (right-hand side) of the
     610             :  * equation which determines the radial (y) dependence of the Bondi quantity
     611             :  * \f$H\f$.
     612             :  *
     613             :  * \details The quantity \f$H \equiv \partial_u J\f$ (evaluated at constant y)
     614             :  * is defined via the Bondi form of the metric:
     615             :  * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2
     616             :  * e^{2 \beta} du dr - 2 r^2 h_{AB} U^B du dx^A + r^2 h_{A B} dx^A dx^B. \f]
     617             :  * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical
     618             :  * angular dyad \f$q^A\f$:
     619             :  * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B.\f]
     620             :  * See \cite Bishop1997ik \cite Handmer2014qha for full details.
     621             :  *
     622             :  * We write the equations of motion in the compactified coordinate \f$ y \equiv
     623             :  * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the
     624             :  * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of
     625             :  * the worldtube. The equation which determines \f$W\f$ on a surface of constant
     626             :  * \f$u\f$ given \f$J\f$,\f$\beta\f$, \f$Q\f$, \f$U\f$, and \f$W\f$ on the same
     627             :  * surface is written as
     628             :  * \f[(1 - y) \partial_y H + H + (1 - y)(\mathcal{D}_J H
     629             :  * + \bar{\mathcal{D}}_J \bar{H}) = A_J + (1 - y) B_J.\f]
     630             :  *
     631             :  * We refer to \f$A_J\f$ as the "pole part" of the integrand
     632             :  * and \f$B_J\f$ as the "regular part". The pole part is computed by this
     633             :  * function, and has the expression
     634             :  * \f{align*}
     635             :  * A_J =& - \tfrac{1}{2} \eth (J \bar{U}) -  \eth (\bar{U}) J -  \tfrac{1}{2}
     636             :  * \bar{\eth} (U) J -  \eth (U) K -  \tfrac{1}{2} (\bar{\eth} (J) U) + 2 J W
     637             :  * \f}
     638             :  */
     639             : template <>
     640           1 : struct ComputeBondiIntegrand<Tags::PoleOfIntegrand<Tags::BondiH>> {
     641             :  public:
     642           0 :   using pre_swsh_derivative_tags =
     643             :       tmpl::list<Tags::BondiJ, Tags::BondiU, Tags::BondiW>;
     644           0 :   using swsh_derivative_tags = tmpl::list<
     645             :       Spectral::Swsh::Tags::Derivative<Tags::BondiU, Spectral::Swsh::Tags::Eth>,
     646             :       Spectral::Swsh::Tags::Derivative<Tags::BondiJ,
     647             :                                        Spectral::Swsh::Tags::Ethbar>,
     648             :       Spectral::Swsh::Tags::Derivative<
     649             :           ::Tags::Multiplies<Tags::BondiJbar, Tags::BondiU>,
     650             :           Spectral::Swsh::Tags::Ethbar>,
     651             :       Spectral::Swsh::Tags::Derivative<Tags::BondiU,
     652             :                                        Spectral::Swsh::Tags::Ethbar>>;
     653           0 :   using integration_independent_tags = tmpl::list<Tags::BondiK>;
     654           0 :   using temporary_tags = tmpl::list<>;
     655             : 
     656           0 :   using return_tags =
     657             :       tmpl::append<tmpl::list<Tags::PoleOfIntegrand<Tags::BondiH>>,
     658             :                    temporary_tags>;
     659           0 :   using argument_tags =
     660             :       tmpl::append<pre_swsh_derivative_tags, swsh_derivative_tags,
     661             :                    integration_independent_tags>;
     662             : 
     663             :   template <typename... Args>
     664           0 :   static void apply(
     665             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*>
     666             :           pole_of_integrand_for_h,
     667             :       const Args&... args) {
     668             :     apply_impl(make_not_null(&get(*pole_of_integrand_for_h)), get(args)...);
     669             :   }
     670             : 
     671             :  private:
     672           0 :   static void apply_impl(
     673             :       gsl::not_null<SpinWeighted<ComplexDataVector, 2>*>
     674             :           pole_of_integrand_for_h,
     675             :       const SpinWeighted<ComplexDataVector, 2>& j,
     676             :       const SpinWeighted<ComplexDataVector, 1>& u,
     677             :       const SpinWeighted<ComplexDataVector, 0>& w,
     678             :       const SpinWeighted<ComplexDataVector, 2>& eth_u,
     679             :       const SpinWeighted<ComplexDataVector, 1>& ethbar_j,
     680             :       const SpinWeighted<ComplexDataVector, -2>& ethbar_jbar_u,
     681             :       const SpinWeighted<ComplexDataVector, 0>& ethbar_u,
     682             :       const SpinWeighted<ComplexDataVector, 0>& k);
     683             : };
     684             : 
     685             : /*!
     686             :  * \brief Computes the pole part of the integrand (right-hand side) of the
     687             :  * equation which determines the radial (y) dependence of the Bondi quantity
     688             :  * \f$H\f$.
     689             :  *
     690             :  * \details The quantity \f$H \equiv \partial_u J\f$ (evaluated at constant y)
     691             :  * is defined via the Bondi form of the metric:
     692             :  * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2
     693             :  * e^{2 \beta} du dr - 2 r^2 h_{AB} U^B du dx^A + r^2 h_{A B} dx^A dx^B. \f]
     694             :  * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical
     695             :  * angular dyad \f$q^A\f$:
     696             :  * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B.\f]
     697             :  * See \cite Bishop1997ik \cite Handmer2014qha for full details.
     698             :  *
     699             :  * We write the equations of motion in the compactified coordinate \f$ y \equiv
     700             :  * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the
     701             :  * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of
     702             :  * the worldtube. The equation which determines \f$H\f$ on a surface of constant
     703             :  * \f$u\f$ given \f$J\f$,\f$\beta\f$, \f$Q\f$, \f$U\f$, and \f$W\f$ on the same
     704             :  * surface is written as
     705             :  * \f[(1 - y) \partial_y H + H + (1 - y)(\mathcal{D}_J H + \bar{\mathcal{D}}_J
     706             :  * \bar{H}) = A_J + (1 - y) B_J.\f]
     707             :  * We refer to \f$A_J\f$ as the "pole part" of the integrand
     708             :  * and \f$B_J\f$ as the "regular part". The pole part is computed by this
     709             :  * function, and has the expression
     710             :  * \f{align*}
     711             :  * B_J =& -\tfrac{1}{2} \left(\eth(\partial_y (J) \bar{U}) +  \partial_y
     712             :  * (\bar{\eth} (J)) U \right) + J (\mathcal{B}_J + \bar{\mathcal{B}}_J) \notag\\
     713             :  * &+ \frac{e^{2 \beta}}{2 R} \left(\mathcal{C}_J + \eth (\eth (\beta)) -
     714             :  * \tfrac{1}{2} \eth (Q) -  (\mathcal{A}_J + \bar{\mathcal{A}_J}) J +
     715             :  * \frac{\bar{\mathcal{C}_J} J^2}{K^2} + \frac{\eth (J (-2 \bar{\eth} (\beta) +
     716             :  * \bar{Q}))}{4 K} -  \frac{\eth (\bar{Q}) J}{4 K} + (\eth (\beta) -
     717             :  * \tfrac{1}{2} Q)^2\right) \notag\\
     718             :  * &-  \partial_y (J)  \left(\frac{\eth (U) \bar{J}}{2 K} -  \tfrac{1}{2}
     719             :  * \bar{\eth} (\bar{U}) J K + \tfrac{1}{4} (\eth (\bar{U}) -  \bar{\eth} (U))
     720             :  * K^2 + \frac{1}{2} \frac{\eth (R) \bar{U}}{R} -  \frac{1}{2}
     721             :  * W\right)\notag\\
     722             :  * &+  \partial_y (\bar{J}) \left(- \tfrac{1}{4} (- \eth (\bar{U}) + \bar{\eth}
     723             :  * (U)) J^2 + \eth (U) J \left(- \frac{1}{2 K} + \tfrac{1}{2}
     724             :  * K\right)\right)\notag\\
     725             :  * &+ (1 - y) \bigg[\frac{1}{2} \left(- \frac{\partial_y (J)}{R} + \frac{2
     726             :  * \partial_{u} (R) \partial_y (\partial_y (J))}{R} + \partial_y
     727             :  * (\partial_y (J)) W\right)  + \partial_y (J) \left(\tfrac{1}{2} \partial_y (W)
     728             :  * + \frac{1}{2 R}\right)\bigg]\notag\\
     729             :  * &+ (1 - y)^2 \bigg[ \frac{\partial_y (\partial_y (J)) }{4 R} \bigg],
     730             :  * \f}
     731             :  * where
     732             :  * \f{align*}
     733             :  * \mathcal{A}_J =& \tfrac{1}{4} \eth (\eth (\bar{J})) -  \frac{1}{4 K^3} -
     734             :  * \frac{\eth (\bar{\eth} (J \bar{J})) -  (\eth (\bar{\eth} (\bar{J})) - 4
     735             :  \bar{J})
     736             :  * J}{16 K^3} + \frac{3}{4 K} -  \frac{\eth (\bar{\eth} (\beta))}{4 K} \notag\\
     737             :  * &-  \frac{\eth (\bar{\eth} (J)) \bar{J} (1 -  \frac{1}{4 K^2})}{4 K} +
     738             :  * \tfrac{1}{2} \eth (\bar{J}) \left(\eth (\beta) + \frac{\bar{\eth} (J \bar{J})
     739             :  * J}{4 K^3} -  \frac{\bar{\eth} (J) (-1 + 2 K^2)}{4 K^3} -  \tfrac{1}{2}
     740             :  * Q\right)\\
     741             :  * \mathcal{B}_J =& - \frac{\eth (U) \bar{J} \partial_y (J \bar{J})}{4 K} +
     742             :  * \tfrac{1}{2} \partial_y (W) + \frac{1}{4 R} + \tfrac{1}{4} \bar{\eth} (J)
     743             :  * \partial_y (\bar{J}) U -  \frac{\bar{\eth} (J \bar{J}) \partial_y (J \bar{J})
     744             :  * U}{8 K^2} \notag\\&-  \tfrac{1}{4} J \partial_y (\eth (\bar{J})) \bar{U} +
     745             :  * \tfrac{1}{4} (\eth (J \partial_y (\bar{J})) + \frac{J \eth (R)
     746             :  * \partial_y(\bar{J})}{R}) \bar{U} \\
     747             :  * &+ (1 - y) \bigg[ \frac{\mathcal{D}_J \partial_{u} (R)
     748             :  * \partial_y (J)}{R} -  \tfrac{1}{4} \partial_y (J) \partial_y (\bar{J}) W +
     749             :  * \frac{(\partial_y (J \bar{J}))^2 W}{16 K^2} \bigg] \\
     750             :  * &+ (1 - y)^2 \bigg[ - \frac{\partial_y (J) \partial_y (\bar{J})}{8 R} +
     751             :  * \frac{(\partial_y (J \bar{J}))^2}{32 K^2 R} \bigg]\\
     752             :  * \mathcal{C}_J =& \tfrac{1}{2} \bar{\eth} (J) K (\eth (\beta) -  \tfrac{1}{2}
     753             :  * Q)\\ \mathcal{D}_J =& \tfrac{1}{4} \left(-2 \partial_y (\bar{J}) +
     754             :  * \frac{\bar{J} \partial_y (J \bar{J})}{K^2}\right)
     755             :  * \f}
     756             :  */
     757             : template <>
     758           1 : struct ComputeBondiIntegrand<Tags::RegularIntegrand<Tags::BondiH>> {
     759             :  public:
     760           0 :   using pre_swsh_derivative_tags =
     761             :       tmpl::list<Tags::Dy<Tags::Dy<Tags::BondiJ>>, Tags::Dy<Tags::BondiJ>,
     762             :                  Tags::Dy<Tags::BondiW>, Tags::Exp2Beta, Tags::BondiJ,
     763             :                  Tags::BondiQ, Tags::BondiU, Tags::BondiW>;
     764           0 :   using swsh_derivative_tags = tmpl::list<
     765             :       Spectral::Swsh::Tags::Derivative<Tags::BondiBeta,
     766             :                                        Spectral::Swsh::Tags::Eth>,
     767             :       Spectral::Swsh::Tags::Derivative<Tags::BondiBeta,
     768             :                                        Spectral::Swsh::Tags::EthEth>,
     769             :       Spectral::Swsh::Tags::Derivative<Tags::BondiBeta,
     770             :                                        Spectral::Swsh::Tags::EthEthbar>,
     771             :       Spectral::Swsh::Tags::Derivative<
     772             :           Spectral::Swsh::Tags::Derivative<Tags::BondiJ,
     773             :                                            Spectral::Swsh::Tags::Ethbar>,
     774             :           Spectral::Swsh::Tags::Eth>,
     775             :       Spectral::Swsh::Tags::Derivative<
     776             :           ::Tags::Multiplies<Tags::BondiJ, Tags::BondiJbar>,
     777             :           Spectral::Swsh::Tags::EthEthbar>,
     778             :       Spectral::Swsh::Tags::Derivative<
     779             :           ::Tags::Multiplies<Tags::BondiJ, Tags::BondiJbar>,
     780             :           Spectral::Swsh::Tags::Eth>,
     781             :       Spectral::Swsh::Tags::Derivative<Tags::BondiQ, Spectral::Swsh::Tags::Eth>,
     782             :       Spectral::Swsh::Tags::Derivative<Tags::BondiU, Spectral::Swsh::Tags::Eth>,
     783             :       Spectral::Swsh::Tags::Derivative<
     784             :           ::Tags::Multiplies<Tags::BondiUbar, Tags::Dy<Tags::BondiJ>>,
     785             :           Spectral::Swsh::Tags::Eth>,
     786             :       Spectral::Swsh::Tags::Derivative<Tags::Dy<Tags::BondiJ>,
     787             :                                        Spectral::Swsh::Tags::Ethbar>,
     788             :       Spectral::Swsh::Tags::Derivative<Tags::BondiJ,
     789             :                                        Spectral::Swsh::Tags::EthbarEthbar>,
     790             :       Spectral::Swsh::Tags::Derivative<Tags::BondiJ,
     791             :                                        Spectral::Swsh::Tags::Ethbar>,
     792             :       Spectral::Swsh::Tags::Derivative<
     793             :           ::Tags::Multiplies<Tags::BondiJbar, Tags::Dy<Tags::BondiJ>>,
     794             :           Spectral::Swsh::Tags::Ethbar>,
     795             :       Spectral::Swsh::Tags::Derivative<Tags::JbarQMinus2EthBeta,
     796             :                                        Spectral::Swsh::Tags::Ethbar>,
     797             :       Spectral::Swsh::Tags::Derivative<Tags::BondiQ,
     798             :                                        Spectral::Swsh::Tags::Ethbar>,
     799             :       Spectral::Swsh::Tags::Derivative<Tags::BondiU,
     800             :                                        Spectral::Swsh::Tags::Ethbar>>;
     801           0 :   using integration_independent_tags =
     802             :       tmpl::list<Tags::DuRDividedByR, Tags::EthRDividedByR, Tags::BondiK,
     803             :                  Tags::OneMinusY, Tags::BondiR>;
     804           0 :   using temporary_tags =
     805             :       tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
     806             :                                       std::integral_constant<int, 0>>,
     807             :                  ::Tags::SpinWeighted<::Tags::TempScalar<1, ComplexDataVector>,
     808             :                                       std::integral_constant<int, 0>>,
     809             :                  ::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
     810             :                                       std::integral_constant<int, 2>>>;
     811             : 
     812           0 :   using return_tags =
     813             :       tmpl::append<tmpl::list<Tags::RegularIntegrand<Tags::BondiH>>,
     814             :                    temporary_tags>;
     815           0 :   using argument_tags =
     816             :       tmpl::append<pre_swsh_derivative_tags, swsh_derivative_tags,
     817             :                    integration_independent_tags>;
     818             : 
     819             :   template <typename... Args>
     820           0 :   static void apply(
     821             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*>
     822             :           regular_integrand_for_h,
     823             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
     824             :           script_aj,
     825             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
     826             :           script_bj,
     827             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*>
     828             :           script_cj,
     829             :       const Args&... args) {
     830             :     apply_impl(make_not_null(&get(*regular_integrand_for_h)),
     831             :                make_not_null(&get(*script_aj)), make_not_null(&get(*script_bj)),
     832             :                make_not_null(&get(*script_cj)), get(args)...);
     833             :   }
     834             : 
     835             :  private:
     836           0 :   static void apply_impl(
     837             :       gsl::not_null<SpinWeighted<ComplexDataVector, 2>*>
     838             :           regular_integrand_for_h,
     839             :       gsl::not_null<SpinWeighted<ComplexDataVector, 0>*> script_aj,
     840             :       gsl::not_null<SpinWeighted<ComplexDataVector, 0>*> script_bj,
     841             :       gsl::not_null<SpinWeighted<ComplexDataVector, 2>*> script_cj,
     842             :       const SpinWeighted<ComplexDataVector, 2>& dy_dy_j,
     843             :       const SpinWeighted<ComplexDataVector, 2>& dy_j,
     844             :       const SpinWeighted<ComplexDataVector, 0>& dy_w,
     845             :       const SpinWeighted<ComplexDataVector, 0>& exp_2_beta,
     846             :       const SpinWeighted<ComplexDataVector, 2>& j,
     847             :       const SpinWeighted<ComplexDataVector, 1>& q,
     848             :       const SpinWeighted<ComplexDataVector, 1>& u,
     849             :       const SpinWeighted<ComplexDataVector, 0>& w,
     850             :       const SpinWeighted<ComplexDataVector, 1>& eth_beta,
     851             :       const SpinWeighted<ComplexDataVector, 2>& eth_eth_beta,
     852             :       const SpinWeighted<ComplexDataVector, 0>& eth_ethbar_beta,
     853             :       const SpinWeighted<ComplexDataVector, 2>& eth_ethbar_j,
     854             :       const SpinWeighted<ComplexDataVector, 0>& eth_ethbar_j_jbar,
     855             :       const SpinWeighted<ComplexDataVector, 1>& eth_j_jbar,
     856             :       const SpinWeighted<ComplexDataVector, 2>& eth_q,
     857             :       const SpinWeighted<ComplexDataVector, 2>& eth_u,
     858             :       const SpinWeighted<ComplexDataVector, 2>& eth_ubar_dy_j,
     859             :       const SpinWeighted<ComplexDataVector, 1>& ethbar_dy_j,
     860             :       const SpinWeighted<ComplexDataVector, 0>& ethbar_ethbar_j,
     861             :       const SpinWeighted<ComplexDataVector, 1>& ethbar_j,
     862             :       const SpinWeighted<ComplexDataVector, -1>& ethbar_jbar_dy_j,
     863             :       const SpinWeighted<ComplexDataVector, -2>& ethbar_jbar_q_minus_2_eth_beta,
     864             :       const SpinWeighted<ComplexDataVector, 0>& ethbar_q,
     865             :       const SpinWeighted<ComplexDataVector, 0>& ethbar_u,
     866             :       const SpinWeighted<ComplexDataVector, 0>& du_r_divided_by_r,
     867             :       const SpinWeighted<ComplexDataVector, 1>& eth_r_divided_by_r,
     868             :       const SpinWeighted<ComplexDataVector, 0>& k,
     869             :       const SpinWeighted<ComplexDataVector, 0>& one_minus_y,
     870             :       const SpinWeighted<ComplexDataVector, 0>& r);
     871             : };
     872             : 
     873             : /*!
     874             :  * \brief Computes the linear factor which multiplies \f$H\f$ in the
     875             :  * equation which determines the radial (y) dependence of the Bondi quantity
     876             :  * \f$H\f$.
     877             :  *
     878             :  * \details The quantity \f$H \equiv \partial_u J\f$ (evaluated at constant y)
     879             :  * is defined via the Bondi form of the metric:
     880             :  * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2
     881             :  * e^{2 \beta} du dr - 2 r^2 h_{AB} U^B du dx^A + r^2 h_{A B} dx^A dx^B. \f]
     882             :  * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical
     883             :  * angular dyad \f$q^A\f$:
     884             :  * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B.\f]
     885             :  * See \cite Bishop1997ik \cite Handmer2014qha for full details.
     886             :  *
     887             :  * We write the equations of motion in the compactified coordinate \f$ y \equiv
     888             :  * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the
     889             :  * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of
     890             :  * the worldtube. The equation which determines \f$H\f$ on a surface of constant
     891             :  * \f$u\f$ given \f$J\f$,\f$\beta\f$, \f$Q\f$, \f$U\f$, and \f$W\f$ on the same
     892             :  * surface is written as
     893             :  * \f[(1 - y) \partial_y H + H + (1 - y) J (\mathcal{D}_J
     894             :  * H + \bar{\mathcal{D}}_J \bar{H}) = A_J + (1 - y) B_J.\f]
     895             :  * The quantity \f$1 +(1 - y) J \mathcal{D}_J\f$ is the linear factor
     896             :  * for the non-conjugated \f$H\f$, and is computed from the equation:
     897             :  * \f[\mathcal{D}_J = \frac{1}{4}(-2 \partial_y \bar{J} + \frac{\bar{J}
     898             :  * \partial_y (J \bar{J})}{K^2})\f]
     899             :  */
     900             : template <>
     901           1 : struct ComputeBondiIntegrand<Tags::LinearFactor<Tags::BondiH>> {
     902             :  public:
     903           0 :   using pre_swsh_derivative_tags =
     904             :       tmpl::list<Tags::Dy<Tags::BondiJ>, Tags::BondiJ>;
     905           0 :   using swsh_derivative_tags = tmpl::list<>;
     906           0 :   using integration_independent_tags = tmpl::list<Tags::OneMinusY>;
     907           0 :   using temporary_tags =
     908             :       tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
     909             :                                       std::integral_constant<int, 2>>>;
     910             : 
     911           0 :   using return_tags = tmpl::append<tmpl::list<Tags::LinearFactor<Tags::BondiH>>,
     912             :                                    temporary_tags>;
     913           0 :   using argument_tags =
     914             :       tmpl::append<pre_swsh_derivative_tags, swsh_derivative_tags,
     915             :                    integration_independent_tags>;
     916             : 
     917             :   template <typename... Args>
     918           0 :   static void apply(
     919             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
     920             :           linear_factor_for_h,
     921             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*>
     922             :           script_djbar,
     923             :       const Args&... args) {
     924             :     apply_impl(make_not_null(&get(*linear_factor_for_h)),
     925             :                make_not_null(&get(*script_djbar)), get(args)...);
     926             :   }
     927             : 
     928             :  private:
     929           0 :   static void apply_impl(
     930             :       gsl::not_null<SpinWeighted<ComplexDataVector, 0>*> linear_factor_for_h,
     931             :       gsl::not_null<SpinWeighted<ComplexDataVector, 2>*> script_djbar,
     932             :       const SpinWeighted<ComplexDataVector, 2>& dy_j,
     933             :       const SpinWeighted<ComplexDataVector, 2>& j,
     934             :       const SpinWeighted<ComplexDataVector, 0>& one_minus_y);
     935             : };
     936             : 
     937             : /*!
     938             :  * \brief Computes the linear factor which multiplies \f$\bar{H}\f$ in the
     939             :  * equation which determines the radial (y) dependence of the Bondi quantity
     940             :  * \f$H\f$.
     941             :  *
     942             :  * \details The quantity \f$H \equiv \partial_u J\f$ (evaluated at constant y)
     943             :  * is defined via the Bondi form of the metric:
     944             :  * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2
     945             :  * e^{2 \beta} du dr - 2 r^2 h_{AB} U^B du dx^A + r^2 h_{A B} dx^A dx^B. \f]
     946             :  * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical
     947             :  * angular dyad \f$q^A\f$:
     948             :  * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B.\f]
     949             :  * See \cite Bishop1997ik \cite Handmer2014qha for full details.
     950             :  *
     951             :  * We write the equations of motion in the compactified coordinate \f$ y \equiv
     952             :  * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the
     953             :  * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of
     954             :  * the worldtube. The equation which determines \f$H\f$ on a surface of constant
     955             :  * \f$u\f$ given \f$J\f$,\f$\beta\f$, \f$Q\f$, \f$U\f$, and \f$W\f$ on the same
     956             :  * surface is written as
     957             :  * \f[(1 - y) \partial_y H + H + (1 - y) J (\mathcal{D}_J H +
     958             :  * \bar{\mathcal{D}}_J \bar{H}) = A_J + (1 - y) B_J.\f]
     959             :  * The quantity \f$ (1 - y) J \bar{\mathcal{D}}_J\f$ is the linear factor
     960             :  * for the non-conjugated \f$H\f$, and is computed from the equation:
     961             :  * \f[\mathcal{D}_J = \frac{1}{4}(-2 \partial_y \bar{J} + \frac{\bar{J}
     962             :  * \partial_y (J \bar{J})}{K^2})\f]
     963             :  */
     964             : template <>
     965           1 : struct ComputeBondiIntegrand<Tags::LinearFactorForConjugate<Tags::BondiH>> {
     966             :  public:
     967           0 :   using pre_swsh_derivative_tags =
     968             :       tmpl::list<Tags::Dy<Tags::BondiJ>, Tags::BondiJ>;
     969           0 :   using swsh_derivative_tags = tmpl::list<>;
     970           0 :   using integration_independent_tags = tmpl::list<Tags::OneMinusY>;
     971           0 :   using temporary_tags =
     972             :       tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>,
     973             :                                       std::integral_constant<int, 2>>>;
     974             : 
     975           0 :   using return_tags =
     976             :       tmpl::append<tmpl::list<Tags::LinearFactorForConjugate<Tags::BondiH>>,
     977             :                    temporary_tags>;
     978           0 :   using argument_tags =
     979             :       tmpl::append<pre_swsh_derivative_tags, swsh_derivative_tags,
     980             :                    integration_independent_tags>;
     981             : 
     982             :   template <typename... Args>
     983           0 :   static void apply(
     984             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 4>>*>
     985             :           linear_factor_for_conjugate_h,
     986             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*>
     987             :           script_djbar,
     988             :       const Args&... args) {
     989             :     apply_impl(make_not_null(&get(*linear_factor_for_conjugate_h)),
     990             :                make_not_null(&get(*script_djbar)), get(args)...);
     991             :   }
     992             : 
     993             :  private:
     994           0 :   static void apply_impl(
     995             :       gsl::not_null<SpinWeighted<ComplexDataVector, 4>*>
     996             :           linear_factor_for_conjugate_h,
     997             :       gsl::not_null<SpinWeighted<ComplexDataVector, 2>*> script_djbar,
     998             :       const SpinWeighted<ComplexDataVector, 2>& dy_j,
     999             :       const SpinWeighted<ComplexDataVector, 2>& j,
    1000             :       const SpinWeighted<ComplexDataVector, 0>& one_minus_y);
    1001             : };
    1002             : 
    1003             : /*!
    1004             :  *\brief Computes the pole part of the integrand (right-hand side) of the
    1005             :  * equation which determines the radial (y) dependence of the scalar quantity
    1006             :  * \f$\Pi\f$.
    1007             :  *
    1008             :  * \details The evolution equation for \f$\Pi\f$ is written as
    1009             :  * \f[(1 - y) \partial_y \Pi + \Pi = A_\Pi + (1 - y) B_\Pi.\f]
    1010             :  * We refer to \f$A_\Pi\f$ as the "pole part" of the integrand and \f$B_\Pi\f$
    1011             :  * as the "regular part". The pole part is computed by this function, and has
    1012             :  * the expression
    1013             :  * \f[A_\Pi = - \Re \left(U \bar{\eth}\psi\right).\f]
    1014             :  */
    1015             : template <>
    1016           1 : struct ComputeBondiIntegrand<Tags::PoleOfIntegrand<Tags::KleinGordonPi>> {
    1017             :  public:
    1018           0 :   using pre_swsh_derivative_tags = tmpl::list<>;
    1019           0 :   using swsh_derivative_tags =
    1020             :       tmpl::list<Spectral::Swsh::Tags::Derivative<Tags::KleinGordonPsi,
    1021             :                                                   Spectral::Swsh::Tags::Eth>>;
    1022           0 :   using integration_independent_tags = tmpl::list<Tags::BondiU>;
    1023             : 
    1024           0 :   using return_tags = tmpl::list<Tags::PoleOfIntegrand<Tags::KleinGordonPi>>;
    1025           0 :   using argument_tags =
    1026             :       tmpl::append<pre_swsh_derivative_tags, swsh_derivative_tags,
    1027             :                    integration_independent_tags>;
    1028             : 
    1029             :   template <typename... Args>
    1030           0 :   static void apply(
    1031             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
    1032             :           pole_of_integrand_for_kg_pi,
    1033             :       const Args&... args) {
    1034             :     apply_impl(make_not_null(&get(*pole_of_integrand_for_kg_pi)), get(args)...);
    1035             :   }
    1036             : 
    1037             :  private:
    1038           0 :   static void apply_impl(gsl::not_null<SpinWeighted<ComplexDataVector, 0>*>
    1039             :                              pole_of_integrand_for_kg_pi,
    1040             :                          const SpinWeighted<ComplexDataVector, 1>& eth_kg_psi,
    1041             :                          const SpinWeighted<ComplexDataVector, 1>& bondi_u);
    1042             : };
    1043             : 
    1044             : /*!
    1045             :  * \brief Computes the regular part of the integrand (right-hand side) of the
    1046             :  * equation which determines the radial (y) dependence of the scalar quantity
    1047             :  * \f$\Pi\f$.
    1048             :  *
    1049             :  * \details The evolution equation for \f$\Pi\f$ is written as
    1050             :  * \f[(1 - y) \partial_y \Pi + \Pi = A_\Pi + (1 - y) B_\Pi.\f]
    1051             :  * We refer to \f$A_\Pi\f$ as the "pole part" of the integrand and \f$B_\Pi\f$
    1052             :  * as the "regular part". The regular part is computed by this function, and has
    1053             :  * the expression:
    1054             :  * \f[ B_\Pi = \frac{1}{4R}e^{2\beta}\left(N_{\psi 1}-N_{\psi 2} +N_{\psi
    1055             :  * 3}\right) -\frac{R}{2}\frac{N_{\psi 4}}{(1-y)^2}
    1056             :  * +\tau+\frac{\partial_{u}R}{R}(1-y)\partial_y^2\psi, \f]
    1057             :  * where
    1058             :  * \f{align*}{
    1059             :  * & N_{\psi 1}= K(2\Re \eth\beta\bar{\eth}\psi + \eth\bar{\eth}\psi) \\
    1060             :  * & N_{\psi 2}=\Re \bar{J}\eth\eth\psi + 2\Re \bar{\eth}\beta \bar{\eth}\psi
    1061             :  * J + \Re \bar{\eth}J\bar{\eth}\psi \\
    1062             :  * & N_{\psi 3}=\Re \eth K\bar{\eth}\psi \\
    1063             :  * & N_{\psi 4}=2\Re (\bar{\eth}\psi ) \partial_rU + 2\Re \eth
    1064             :  * \bar{U}\partial_r\psi + 4\Re \bar{U}\eth\partial_r\psi \\
    1065             :  * & \tau=\frac{1}{2}(1-y)\partial_y W\partial_y\psi +
    1066             :  * \frac{(1-y)^2}{4R}\partial_y^2\psi +
    1067             :  * \frac{1}{2}(1-y)W\partial_y^2\psi + \frac{1}{2}W\partial_y\psi
    1068             :  * \f}
    1069             :  */
    1070             : template <>
    1071           1 : struct ComputeBondiIntegrand<Tags::RegularIntegrand<Tags::KleinGordonPi>> {
    1072             :  public:
    1073           0 :   using pre_swsh_derivative_tags =
    1074             :       tmpl::list<Tags::Dy<Tags::Dy<Tags::KleinGordonPsi>>,
    1075             :                  Tags::Dy<Tags::KleinGordonPsi>, Tags::Dy<Tags::BondiU>,
    1076             :                  Tags::Dy<Tags::BondiW>>;
    1077           0 :   using swsh_derivative_tags =
    1078             :       tmpl::list<Spectral::Swsh::Tags::Derivative<
    1079             :                      Tags::Dy<Tags::KleinGordonPsi>, Spectral::Swsh::Tags::Eth>,
    1080             :                  Spectral::Swsh::Tags::Derivative<Tags::KleinGordonPsi,
    1081             :                                                   Spectral::Swsh::Tags::EthEth>,
    1082             :                  Spectral::Swsh::Tags::Derivative<Tags::KleinGordonPsi,
    1083             :                                                   Spectral::Swsh::Tags::Eth>,
    1084             :                  Spectral::Swsh::Tags::Derivative<Tags::BondiJ,
    1085             :                                                   Spectral::Swsh::Tags::Ethbar>,
    1086             :                  Spectral::Swsh::Tags::Derivative<Tags::BondiU,
    1087             :                                                   Spectral::Swsh::Tags::Ethbar>,
    1088             :                  Spectral::Swsh::Tags::Derivative<Tags::BondiBeta,
    1089             :                                                   Spectral::Swsh::Tags::Eth>,
    1090             :                  Spectral::Swsh::Tags::Derivative<
    1091             :                      ::Tags::Multiplies<Tags::BondiJ, Tags::BondiJbar>,
    1092             :                      Spectral::Swsh::Tags::Eth>,
    1093             :                  Spectral::Swsh::Tags::Derivative<
    1094             :                      Tags::KleinGordonPsi, Spectral::Swsh::Tags::EthEthbar>>;
    1095           0 :   using integration_independent_tags =
    1096             :       tmpl::list<Tags::BondiJ, Tags::Exp2Beta, Tags::DuRDividedByR,
    1097             :                  Tags::OneMinusY, Tags::EthRDividedByR, Tags::BondiR,
    1098             :                  Tags::BondiK, Tags::BondiU, Tags::BondiW>;
    1099             : 
    1100           0 :   using return_tags = tmpl::list<Tags::RegularIntegrand<Tags::KleinGordonPi>>;
    1101           0 :   using argument_tags =
    1102             :       tmpl::append<pre_swsh_derivative_tags, swsh_derivative_tags,
    1103             :                    integration_independent_tags>;
    1104             : 
    1105             :   template <typename... Args>
    1106           0 :   static void apply(
    1107             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*>
    1108             :           regular_integrand_for_kg_pi,
    1109             :       const Args&... args) {
    1110             :     apply_impl(make_not_null(&get(*regular_integrand_for_kg_pi)), get(args)...);
    1111             :   }
    1112             : 
    1113             :  private:
    1114           0 :   static void apply_impl(
    1115             :       gsl::not_null<SpinWeighted<ComplexDataVector, 0>*>
    1116             :           regular_integrand_for_kg_pi,
    1117             :       // pre_swsh_derivative_tags
    1118             :       const SpinWeighted<ComplexDataVector, 0>& dy_dy_kg_psi,
    1119             :       const SpinWeighted<ComplexDataVector, 0>& dy_kg_psi,
    1120             :       const SpinWeighted<ComplexDataVector, 1>& dy_bondi_u,
    1121             :       const SpinWeighted<ComplexDataVector, 0>& dy_w,
    1122             :       // swsh_derivative_tags
    1123             :       const SpinWeighted<ComplexDataVector, 1>& eth_dy_kg_psi,
    1124             :       const SpinWeighted<ComplexDataVector, 2>& eth_eth_kg_psi,
    1125             :       const SpinWeighted<ComplexDataVector, 1>& eth_kg_psi,
    1126             :       const SpinWeighted<ComplexDataVector, 1>& ethbar_j,
    1127             :       const SpinWeighted<ComplexDataVector, 0>& ethbar_u,
    1128             :       const SpinWeighted<ComplexDataVector, 1>& eth_beta,
    1129             :       const SpinWeighted<ComplexDataVector, 1>& eth_j_jbar,
    1130             :       const SpinWeighted<ComplexDataVector, 0>& eth_ethbar_kg_psi,
    1131             :       // integration_independent_tags
    1132             :       const SpinWeighted<ComplexDataVector, 2>& j,
    1133             :       const SpinWeighted<ComplexDataVector, 0>& exp2beta,
    1134             :       const SpinWeighted<ComplexDataVector, 0>& du_r_divided_by_r,
    1135             :       const SpinWeighted<ComplexDataVector, 0>& one_minus_y,
    1136             :       const SpinWeighted<ComplexDataVector, 1>& eth_r_divided_by_r,
    1137             :       const SpinWeighted<ComplexDataVector, 0>& bondi_r,
    1138             :       const SpinWeighted<ComplexDataVector, 0>& bondi_k,
    1139             :       const SpinWeighted<ComplexDataVector, 1>& bondi_u,
    1140             :       const SpinWeighted<ComplexDataVector, 0>& bondi_w);
    1141             : };
    1142             : }  // namespace Cce

Generated by: LCOV version 1.14