SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/Cce - SwshDerivatives.hpp Hit Total Coverage
Commit: 817e13c5144619b701c7cd870655d8dbf94ab8ce Lines: 9 58 15.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 <cmath>
       7             : #include <cstddef>
       8             : #include <type_traits>
       9             : 
      10             : #include "DataStructures/ComplexDataVector.hpp"
      11             : #include "DataStructures/DataBox/DataBox.hpp"
      12             : #include "DataStructures/Tags.hpp"
      13             : #include "DataStructures/Tensor/Tensor.hpp"
      14             : #include "DataStructures/Tensor/TypeAliases.hpp"
      15             : #include "Evolution/Systems/Cce/IntegrandInputSteps.hpp"
      16             : #include "Evolution/Systems/Cce/Tags.hpp"
      17             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      18             : #include "NumericalAlgorithms/Spectral/Spectral.hpp"
      19             : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshDerivatives.hpp"
      20             : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshTags.hpp"
      21             : #include "Utilities/ForceInline.hpp"
      22             : #include "Utilities/Gsl.hpp"
      23             : #include "Utilities/TypeTraits.hpp"
      24             : #include "Utilities/TypeTraits/IsA.hpp"
      25             : 
      26             : namespace Cce {
      27             : namespace detail {
      28             : // Precomputation routines for supplying the additional quantities necessary
      29             : // to correct the output of the angular derivative routines from the angular
      30             : // derivatives evaluated at constant numerical coordinates (which is what is
      31             : // returned after the libsharp evaluation) to the angular derivatives at
      32             : // constant Bondi radius (which is what appears in the literature equations
      33             : // and is simple to combine to obtain the hypersurface integrands).
      34             : //
      35             : // Warning: this 'on demand' template is a way of taking advantage of the blaze
      36             : // expression templates in a generic, modular way. However, this can be
      37             : // dangerous. The returned value MUST be a blaze expression template directly,
      38             : // and not a wrapper type (like `SpinWeighted`). Otherwise, some information is
      39             : // lost on the stack and the expression template is corrupted. So, in these 'on
      40             : // demand' returns, the arguments must be fully unpacked to vector types.
      41             : //
      42             : // The `Select` operand is a `std::bool_constant` that ensures mutual
      43             : // exclusivity of the template specializations.
      44             : template <typename Tag, typename SpinConstant, typename Select>
      45             : struct OnDemandInputsForSwshJacobianImpl;
      46             : 
      47             : // default to just retrieving it from the box if not providing an expression
      48             : // template shortcut
      49             : template <typename Tag>
      50             : struct OnDemandInputsForSwshJacobianImpl<
      51             :     Tags::Dy<Tag>, std::integral_constant<int, Tag::type::type::spin>,
      52             :     std::bool_constant<not tt::is_a_v<::Tags::Multiplies, Tag> and
      53             :                        not tt::is_a_v<Tags::Dy, Tag>>> {
      54             :   template <typename DataBoxTagList>
      55             :   SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
      56             :       const db::DataBox<DataBoxTagList>& box) {
      57             :     return get(db::get<Tags::Dy<Tag>>(box)).data();
      58             :   }
      59             : };
      60             : 
      61             : // default to retrieving from the box if the requested tag is a second
      62             : // derivative without an additional evaluation channel
      63             : template <typename Tag>
      64             : struct OnDemandInputsForSwshJacobianImpl<
      65             :     Tags::Dy<Tags::Dy<Tag>>, std::integral_constant<int, Tag::type::type::spin>,
      66             :     std::bool_constant<not tt::is_a_v<::Tags::Multiplies, Tag>>> {
      67             :   template <typename DataBoxTagList>
      68             :   SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
      69             :       const db::DataBox<DataBoxTagList>& box) {
      70             :     return get(db::get<Tags::Dy<Tags::Dy<Tag>>>(box)).data();
      71             :   }
      72             : };
      73             : 
      74             : // use the product rule to provide an expression template for derivatives of
      75             : // products
      76             : template <typename LhsTag, typename RhsTag>
      77             : struct OnDemandInputsForSwshJacobianImpl<
      78             :     Tags::Dy<::Tags::Multiplies<LhsTag, RhsTag>>,
      79             :     std::integral_constant<int,
      80             :                            LhsTag::type::type::spin + RhsTag::type::type::spin>,
      81             :     std::bool_constant<not std::is_same_v<LhsTag, Tags::BondiJbar> and
      82             :                        not std::is_same_v<LhsTag, Tags::BondiUbar> and
      83             :                        not std::is_same_v<RhsTag, Tags::BondiJbar>>> {
      84             :   template <typename DataBoxTagList>
      85             :   SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
      86             :       const db::DataBox<DataBoxTagList>& box) {
      87             :     decltype(auto) lhs = get(db::get<LhsTag>(box)).data();
      88             :     decltype(auto) dy_lhs = get(db::get<Tags::Dy<LhsTag>>(box)).data();
      89             :     decltype(auto) rhs = get(db::get<RhsTag>(box)).data();
      90             :     decltype(auto) dy_rhs = get(db::get<Tags::Dy<RhsTag>>(box)).data();
      91             :     return lhs * dy_rhs + dy_lhs * rhs;
      92             :   }
      93             : };
      94             : 
      95             : // use the product rule and an explicit conjugate to provide an expression
      96             : // template for derivatives of products with `Tags::BondiJbar` as the right-hand
      97             : // operand
      98             : template <typename LhsTag>
      99             : struct OnDemandInputsForSwshJacobianImpl<
     100             :     Tags::Dy<::Tags::Multiplies<LhsTag, Tags::BondiJbar>>,
     101             :     std::integral_constant<int, LhsTag::type::type::spin - 2>, std::true_type> {
     102             :   template <typename DataBoxTagList>
     103             :   SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
     104             :       const db::DataBox<DataBoxTagList>& box) {
     105             :     decltype(auto) lhs = get(get<LhsTag>(box)).data();
     106             :     decltype(auto) dy_lhs = get(get<Tags::Dy<LhsTag>>(box)).data();
     107             :     decltype(auto) jbar = conj(get(get<Tags::BondiJ>(box)).data());
     108             :     decltype(auto) dy_jbar = conj(get(get<Tags::Dy<Tags::BondiJ>>(box)).data());
     109             :     return lhs * dy_jbar + dy_lhs * jbar;
     110             :   }
     111             : };
     112             : 
     113             : // use the product rule and an explicit conjugate to provide an expression
     114             : // template for derivatives of products with `Tags::BondiJbar` as the left-hand
     115             : // operand.
     116             : template <typename RhsTag>
     117             : struct OnDemandInputsForSwshJacobianImpl<
     118             :     Tags::Dy<::Tags::Multiplies<Tags::BondiJbar, RhsTag>>,
     119             :     std::integral_constant<int, RhsTag::type::type::spin - 2>, std::true_type> {
     120             :   template <typename DataBoxTagList>
     121             :   SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
     122             :       const db::DataBox<DataBoxTagList>& box) {
     123             :     decltype(auto) rhs = get(get<RhsTag>(box)).data();
     124             :     decltype(auto) dy_rhs = get(get<Tags::Dy<RhsTag>>(box)).data();
     125             :     decltype(auto) jbar = conj(get(get<Tags::BondiJ>(box)).data());
     126             :     decltype(auto) dy_jbar = conj(get(get<Tags::Dy<Tags::BondiJ>>(box)).data());
     127             :     return dy_jbar * rhs + jbar * dy_rhs;
     128             :   }
     129             : };
     130             : 
     131             : // use the product rule and an explicit conjugate to provide an expression
     132             : // template for derivatives of products with `Tags::BondiUbar` as the left-hand
     133             : // operand.
     134             : template <typename RhsTag>
     135             : struct OnDemandInputsForSwshJacobianImpl<
     136             :     Tags::Dy<::Tags::Multiplies<Tags::BondiUbar, RhsTag>>,
     137             :     std::integral_constant<int, RhsTag::type::type::spin - 1>, std::true_type> {
     138             :   template <typename DataBoxTagList>
     139             :   SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
     140             :       const db::DataBox<DataBoxTagList>& box) {
     141             :     decltype(auto) ubar = conj(get(get<Tags::BondiU>(box)).data());
     142             :     decltype(auto) dy_ubar = conj(get(get<Tags::Dy<Tags::BondiU>>(box)).data());
     143             :     decltype(auto) rhs = get(get<RhsTag>(box)).data();
     144             :     decltype(auto) dy_rhs = get(get<Tags::Dy<RhsTag>>(box)).data();
     145             :     return ubar * dy_rhs + dy_ubar * rhs;
     146             :   }
     147             : };
     148             : 
     149             : // use the product rule and an explicit conjugate to provide an expression
     150             : // template for second derivatives of products with `Tags::BondiJbar` as the
     151             : // right-hand operand.
     152             : template <typename LhsTag>
     153             : struct OnDemandInputsForSwshJacobianImpl<
     154             :     Tags::Dy<Tags::Dy<::Tags::Multiplies<LhsTag, Tags::BondiJbar>>>,
     155             :     std::integral_constant<int, LhsTag::type::type::spin - 2>, std::true_type> {
     156             :   template <typename DataBoxTagList>
     157             :   SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
     158             :       const db::DataBox<DataBoxTagList>& box) {
     159             :     decltype(auto) lhs = get(get<LhsTag>(box)).data();
     160             :     decltype(auto) dy_lhs = get(get<Tags::Dy<LhsTag>>(box)).data();
     161             :     decltype(auto) dy_dy_lhs = get(get<Tags::Dy<Tags::Dy<LhsTag>>>(box)).data();
     162             :     decltype(auto) jbar = conj(get(get<Tags::BondiJ>(box)).data());
     163             :     decltype(auto) dy_jbar = conj(get(get<Tags::Dy<Tags::BondiJ>>(box)).data());
     164             :     decltype(auto) dy_dy_jbar =
     165             :         conj(get(get<Tags::Dy<Tags::Dy<Tags::BondiJ>>>(box)).data());
     166             :     return lhs * dy_dy_jbar + 2.0 * dy_lhs * dy_jbar + dy_dy_lhs * jbar;
     167             :   }
     168             : };
     169             : 
     170             : // default to extracting directly from the box for spin-weighted derivatives of
     171             : // radial (y) derivatives
     172             : template <typename Tag, typename DerivKind>
     173             : struct OnDemandInputsForSwshJacobianImpl<
     174             :     Spectral::Swsh::Tags::Derivative<Tags::Dy<Tag>, DerivKind>,
     175             :     std::integral_constant<
     176             :         int, Spectral::Swsh::Tags::Derivative<Tags::Dy<Tag>, DerivKind>::spin>,
     177             :     std::true_type> {
     178             :   template <typename DataBoxTagList>
     179             :   SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
     180             :       const db::DataBox<DataBoxTagList>& box) {
     181             :     return get(get<Spectral::Swsh::Tags::Derivative<Tags::Dy<Tag>, DerivKind>>(
     182             :                    box))
     183             :         .data();
     184             :   }
     185             : };
     186             : 
     187             : // compute the derivative of the `jbar * (q - 2 eth_beta)` using the product
     188             : // rule and the commutation rule for the partials with respect to y and the
     189             : // spin_weighted derivatives
     190             : template <>
     191             : struct OnDemandInputsForSwshJacobianImpl<Tags::Dy<Tags::JbarQMinus2EthBeta>,
     192             :                                          std::integral_constant<int, -1>,
     193             :                                          std::true_type> {
     194             :   template <typename DataBoxTagList>
     195             :   SPECTRE_ALWAYS_INLINE decltype(auto) operator()(
     196             :       const db::DataBox<DataBoxTagList>& box) {
     197             :     decltype(auto) dy_beta = get(get<Tags::Dy<Tags::BondiBeta>>(box)).data();
     198             :     decltype(auto) dy_j = get(get<Tags::Dy<Tags::BondiJ>>(box)).data();
     199             :     decltype(auto) dy_q = get(get<Tags::Dy<Tags::BondiQ>>(box)).data();
     200             :     decltype(auto) eth_beta =
     201             :         get(get<Spectral::Swsh::Tags::Derivative<Tags::BondiBeta,
     202             :                                                  Spectral::Swsh::Tags::Eth>>(
     203             :                 box))
     204             :             .data();
     205             :     decltype(auto) eth_dy_beta =
     206             :         get(get<Spectral::Swsh::Tags::Derivative<Tags::Dy<Tags::BondiBeta>,
     207             :                                                  Spectral::Swsh::Tags::Eth>>(
     208             :                 box))
     209             :             .data();
     210             :     decltype(auto) eth_r_divided_by_r =
     211             :         get(get<Tags::EthRDividedByR>(box)).data();
     212             :     decltype(auto) j = get(get<Tags::BondiJ>(box)).data();
     213             :     decltype(auto) q = get(get<Tags::BondiQ>(box)).data();
     214             :     return conj(j) * dy_q + conj(dy_j) * q - 2.0 * conj(j) * eth_dy_beta -
     215             :            2.0 * eth_beta * conj(dy_j) -
     216             :            2.0 * conj(j) * eth_r_divided_by_r * dy_beta;
     217             :   }
     218             : };
     219             : }  // namespace detail
     220             : 
     221             : /// Provide an expression template or reference to `Tag`, intended for
     222             : /// situations for which a repeated computation is more
     223             : /// desirable than storing a value in the \ref DataBoxGroup (e.g. for
     224             : /// conjugation and simple product rule expansion).
     225             : template <typename Tag>
     226           1 : using OnDemandInputsForSwshJacobian = detail::OnDemandInputsForSwshJacobianImpl<
     227             :     Tag, std::integral_constant<int, Tag::type::type::spin>, std::true_type>;
     228             : 
     229             : /*!
     230             :  * \brief Performs a mutation to a spin-weighted spherical harmonic derivative
     231             :  * value from the numerical coordinate (the spin-weighted derivative at
     232             :  * fixed \f$y\f$) to the Bondi coordinates (the spin-weighted derivative at
     233             :  * fixed \f$r\f$), inplace to the requested tag.
     234             :  *
     235             :  * \details This should be performed only once for each derivative evaluation
     236             :  * for each tag, as a repeated inplace evaluation will compound and result in
     237             :  * incorrect values in the \ref DataBoxGroup. This is compatible with acting as
     238             :  * a mutation in `db::mutate_apply`.
     239             :  * \note In each specialization, there is an additional type alias
     240             :  * `on_demand_argument_tags` that contains tags that represent additional
     241             :  * quantities to be passed as arguments that need not be in the \ref
     242             :  * DataBoxGroup. These quantities are suggested to be evaluated by the 'on
     243             :  * demand' mechanism provided by `Cce::OnDemandInputsForSwshJacobian`, which
     244             :  * provides the additional quantities as blaze expression templates rather than
     245             :  * unnecessarily caching intermediate results that aren't re-used.
     246             :  */
     247             : template <typename DerivativeTag>
     248           1 : struct ApplySwshJacobianInplace;
     249             : 
     250             : /*!
     251             :  * \brief Specialization for the spin-weighted derivative \f$\eth\f$.
     252             :  *
     253             :  * \details The implemented equation is:
     254             :  *
     255             :  * \f[ \eth F = \eth^\prime F - (1 - y) \frac{\eth R}{R} \partial_y F,
     256             :  * \f]
     257             :  *
     258             :  * where \f$\eth\f$ is the derivative at constant Bondi radius \f$r\f$ and
     259             :  * \f$\eth^\prime\f$ is the derivative at constant numerical radius \f$y\f$.
     260             :  */
     261             : template <typename ArgumentTag>
     262           1 : struct ApplySwshJacobianInplace<
     263             :     Spectral::Swsh::Tags::Derivative<ArgumentTag, Spectral::Swsh::Tags::Eth>> {
     264           0 :   using pre_swsh_derivative_tags = tmpl::list<>;
     265           0 :   using swsh_derivative_tags = tmpl::list<>;
     266           0 :   using integration_independent_tags =
     267             :       tmpl::list<Tags::OneMinusY, Tags::EthRDividedByR>;
     268             : 
     269           0 :   using return_tags = tmpl::list<
     270             :       Spectral::Swsh::Tags::Derivative<ArgumentTag, Spectral::Swsh::Tags::Eth>>;
     271           0 :   using argument_tags = tmpl::append<integration_independent_tags>;
     272           0 :   using on_demand_argument_tags = tmpl::list<Tags::Dy<ArgumentTag>>;
     273             : 
     274           0 :   static constexpr int spin =
     275             :       Spectral::Swsh::Tags::Derivative<ArgumentTag,
     276             :                                        Spectral::Swsh::Tags::Eth>::spin;
     277             :   template <typename DyArgumentType>
     278           0 :   static void apply(
     279             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, spin>>*>
     280             :           eth_argument,
     281             :       const Scalar<SpinWeighted<ComplexDataVector, 0>>& one_minus_y,
     282             :       const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_r_divided_by_r,
     283             :       const DyArgumentType& dy_argument) {
     284             :     get(*eth_argument) -= get(one_minus_y) * get(eth_r_divided_by_r) *
     285             :                           SpinWeighted<DyArgumentType, spin - 1>{dy_argument};
     286             :   }
     287             : };
     288             : 
     289             : /*!
     290             :  * \brief Specialization for the spin-weighted derivative \f$\bar{\eth}\f$.
     291             :  *
     292             :  * \details The implemented equation is:
     293             :  *
     294             :  * \f[
     295             :  * \bar{\eth} F = \bar{\eth}^\prime F
     296             :  * - (1 - y) \frac{\bar{\eth} R}{R} \partial_y F,
     297             :  *\f]
     298             :  *
     299             :  * where \f$\bar{\eth}\f$ is the derivative at constant Bondi radius \f$r\f$ and
     300             :  * \f$\bar{\eth}^\prime\f$ is the derivative at constant numerical radius
     301             :  * \f$y\f$.
     302             :  */
     303             : template <typename ArgumentTag>
     304           1 : struct ApplySwshJacobianInplace<Spectral::Swsh::Tags::Derivative<
     305             :     ArgumentTag, Spectral::Swsh::Tags::Ethbar>> {
     306           0 :   using pre_swsh_derivative_tags = tmpl::list<>;
     307           0 :   using swsh_derivative_tags = tmpl::list<>;
     308           0 :   using integration_independent_tags =
     309             :       tmpl::list<Tags::OneMinusY, Tags::EthRDividedByR>;
     310             : 
     311           0 :   using return_tags = tmpl::list<Spectral::Swsh::Tags::Derivative<
     312             :       ArgumentTag, Spectral::Swsh::Tags::Ethbar>>;
     313           0 :   using argument_tags = tmpl::append<integration_independent_tags>;
     314           0 :   using on_demand_argument_tags = tmpl::list<Tags::Dy<ArgumentTag>>;
     315             : 
     316           0 :   static constexpr int spin =
     317             :       Spectral::Swsh::Tags::Derivative<ArgumentTag,
     318             :                                        Spectral::Swsh::Tags::Ethbar>::spin;
     319             :   template <typename DyArgumentType>
     320           0 :   static void apply(
     321             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, spin>>*>
     322             :           ethbar_argument,
     323             :       const Scalar<SpinWeighted<ComplexDataVector, 0>>& one_minus_y,
     324             :       const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_r_divided_by_r,
     325             :       const DyArgumentType& dy_argument) {
     326             :     get(*ethbar_argument) -=
     327             :         get(one_minus_y) * conj(get(eth_r_divided_by_r)) *
     328             :         SpinWeighted<DyArgumentType, spin + 1>{dy_argument};
     329             :   }
     330             : };
     331             : 
     332             : /*!
     333             :  * \brief Specialization for the spin-weighted derivative \f$\eth \bar{\eth}\f$.
     334             :  *
     335             :  * \details The implemented equation is:
     336             :  *
     337             :  * \f[
     338             :  * \eth \bar{\eth} F = \eth^\prime \bar{\eth}^\prime F
     339             :  * - \frac{\eth R \bar{\eth} R}{R^2} (1 - y)^2 \partial_y^2 F
     340             :  * - (1 - y)\left(\frac{\eth R}{R} \bar{\eth} \partial_y F
     341             :  * + \frac{\bar{\eth} R}{R} \eth \partial_y F
     342             :  * + \frac{\eth \bar\eth R}{R} \partial_y F\right),
     343             :  * \f]
     344             :  *
     345             :  * where \f$\eth \bar{\eth}\f$ is the derivative at constant Bondi radius
     346             :  * \f$r\f$ and \f$\eth^\prime \bar{\eth}^\prime\f$ is the derivative at constant
     347             :  * numerical radius \f$y\f$.
     348             :  */
     349             : template <typename ArgumentTag>
     350           1 : struct ApplySwshJacobianInplace<Spectral::Swsh::Tags::Derivative<
     351             :     ArgumentTag, Spectral::Swsh::Tags::EthEthbar>> {
     352           0 :   using pre_swsh_derivative_tags = tmpl::list<>;
     353           0 :   using swsh_derivative_tags = tmpl::list<>;
     354           0 :   using integration_independent_tags =
     355             :       tmpl::list<Tags::OneMinusY, Tags::EthRDividedByR,
     356             :                  Tags::EthEthbarRDividedByR>;
     357             : 
     358           0 :   using return_tags = tmpl::list<Spectral::Swsh::Tags::Derivative<
     359             :       ArgumentTag, Spectral::Swsh::Tags::EthEthbar>>;
     360           0 :   using argument_tags = tmpl::append<integration_independent_tags>;
     361           0 :   using on_demand_argument_tags =
     362             :       tmpl::list<Tags::Dy<ArgumentTag>, Tags::Dy<Tags::Dy<ArgumentTag>>,
     363             :                  Spectral::Swsh::Tags::Derivative<Tags::Dy<ArgumentTag>,
     364             :                                                   Spectral::Swsh::Tags::Eth>,
     365             :                  Spectral::Swsh::Tags::Derivative<
     366             :                      Tags::Dy<ArgumentTag>, Spectral::Swsh::Tags::Ethbar>>;
     367             : 
     368           0 :   static constexpr int spin =
     369             :       Spectral::Swsh::Tags::Derivative<ArgumentTag,
     370             :                                        Spectral::Swsh::Tags::EthEthbar>::spin;
     371             :   template <typename DyArgumentType, typename DyDyArgumentType,
     372             :             typename EthDyArgumentType, typename EthbarDyArgumentType>
     373           0 :   static void apply(
     374             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, spin>>*>
     375             :           eth_ethbar_argument,
     376             :       const Scalar<SpinWeighted<ComplexDataVector, 0>>& one_minus_y,
     377             :       const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_r_divided_by_r,
     378             :       const Scalar<SpinWeighted<ComplexDataVector, 0>>&
     379             :           eth_ethbar_r_divided_by_r,
     380             :       const DyArgumentType& dy_argument, const DyDyArgumentType& dy_dy_argument,
     381             :       const EthDyArgumentType& eth_dy_argument,
     382             :       const EthbarDyArgumentType ethbar_dy_argument) {
     383             :     get(*eth_ethbar_argument) -=
     384             :         get(eth_r_divided_by_r) * conj(get(eth_r_divided_by_r)) *
     385             :             (square(get(one_minus_y)) *
     386             :              SpinWeighted<DyDyArgumentType, spin>{dy_dy_argument}) +
     387             :         get(one_minus_y) *
     388             :             (get(eth_r_divided_by_r) *
     389             :                  SpinWeighted<EthbarDyArgumentType, spin - 1>{
     390             :                      ethbar_dy_argument} +
     391             :              conj(get(eth_r_divided_by_r)) *
     392             :                  SpinWeighted<EthDyArgumentType, spin + 1>{eth_dy_argument} +
     393             :              get(eth_ethbar_r_divided_by_r) *
     394             :                  SpinWeighted<DyArgumentType, spin>{dy_argument});
     395             :   }
     396             : };
     397             : 
     398             : /*!
     399             :  * \brief Specialization for the spin-weighted derivative \f$\bar{\eth} \eth\f$.
     400             :  *
     401             :  * \details The implemented equation is:
     402             :  *
     403             :  * \f[
     404             :  * \bar{\eth} \eth F = \bar{\eth}^\prime \eth^\prime F
     405             :  * - \frac{\eth R \bar{\eth} R}{R^2} (1 - y)^2 \partial_y^2 F
     406             :  * - (1 - y)\left(\frac{\eth R}{R} \bar{\eth} \partial_y F
     407             :  * + \frac{\bar{\eth} R}{R} \eth \partial_y F
     408             :  * + \frac{\eth \bar\eth R}{R} \partial_y F\right),
     409             :  * \f]
     410             :  *
     411             :  * where \f$\bar{\eth} \eth\f$ is the derivative at constant Bondi radius
     412             :  * \f$r\f$ and \f$\bar{\eth}^\prime \eth^\prime\f$ is the derivative at constant
     413             :  * numerical radius \f$y\f$.
     414             :  */
     415             : template <typename ArgumentTag>
     416           1 : struct ApplySwshJacobianInplace<Spectral::Swsh::Tags::Derivative<
     417             :     ArgumentTag, Spectral::Swsh::Tags::EthbarEth>> {
     418           0 :   using pre_swsh_derivative_tags = tmpl::list<>;
     419           0 :   using swsh_derivative_tags = tmpl::list<>;
     420           0 :   using integration_independent_tags =
     421             :       tmpl::list<Tags::OneMinusY, Tags::EthRDividedByR,
     422             :                  Tags::EthEthbarRDividedByR>;
     423             : 
     424           0 :   using return_tags = tmpl::list<Spectral::Swsh::Tags::Derivative<
     425             :       ArgumentTag, Spectral::Swsh::Tags::EthbarEth>>;
     426           0 :   using argument_tags = tmpl::append<integration_independent_tags>;
     427           0 :   using on_demand_argument_tags =
     428             :       tmpl::list<Tags::Dy<ArgumentTag>, Tags::Dy<Tags::Dy<ArgumentTag>>,
     429             :                  Spectral::Swsh::Tags::Derivative<Tags::Dy<ArgumentTag>,
     430             :                                                   Spectral::Swsh::Tags::Eth>,
     431             :                  Spectral::Swsh::Tags::Derivative<
     432             :                      Tags::Dy<ArgumentTag>, Spectral::Swsh::Tags::Ethbar>>;
     433             : 
     434           0 :   static constexpr int spin =
     435             :       Spectral::Swsh::Tags::Derivative<ArgumentTag,
     436             :                                        Spectral::Swsh::Tags::EthbarEth>::spin;
     437             :   template <typename DyArgumentType, typename DyDyArgumentType,
     438             :             typename EthDyArgumentType, typename EthbarDyArgumentType>
     439           0 :   static void apply(
     440             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, spin>>*>
     441             :           ethbar_eth_argument,
     442             :       const Scalar<SpinWeighted<ComplexDataVector, 0>>& one_minus_y,
     443             :       const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_r_divided_by_r,
     444             :       const Scalar<SpinWeighted<ComplexDataVector, 0>>&
     445             :           eth_ethbar_r_divided_by_r,
     446             :       const DyArgumentType& dy_argument, const DyDyArgumentType& dy_dy_argument,
     447             :       const EthDyArgumentType& eth_dy_argument,
     448             :       const EthbarDyArgumentType ethbar_dy_argument) {
     449             :     get(*ethbar_eth_argument) -=
     450             :         get(eth_r_divided_by_r) * conj(get(eth_r_divided_by_r)) *
     451             :             (square(get(one_minus_y)) *
     452             :              SpinWeighted<DyDyArgumentType, spin>{dy_dy_argument}) +
     453             :         get(one_minus_y) *
     454             :             (get(eth_r_divided_by_r) *
     455             :                  SpinWeighted<EthbarDyArgumentType, spin - 1>{
     456             :                      ethbar_dy_argument} +
     457             :              conj(get(eth_r_divided_by_r)) *
     458             :                  SpinWeighted<EthDyArgumentType, spin + 1>{eth_dy_argument} +
     459             :              get(eth_ethbar_r_divided_by_r) *
     460             :                  SpinWeighted<DyArgumentType, spin>{dy_argument});
     461             :   }
     462             : };
     463             : 
     464             : /*!
     465             :  * \brief Specialization for the spin-weighted derivative \f$\eth \eth\f$.
     466             :  *
     467             :  * \details The implemented equation is:
     468             :  *
     469             :  * \f[
     470             :  * \eth \eth F = \eth^\prime \eth^\prime F
     471             :  * - (1 - y)^2 \frac{(\eth R)^2}{R^2} \partial_y^2 F
     472             :  * - (1 - y) \left( 2 \frac{\eth R}{R} \eth \partial_y F
     473             :  * + \frac{\eth \eth R}{R} \partial_y F\right),
     474             :  * \f]
     475             :  *
     476             :  * where \f$\eth \eth\f$ is the derivative at constant Bondi radius \f$r\f$ and
     477             :  * \f$\eth^\prime \eth^\prime\f$ is the derivative at constant numerical radius
     478             :  * \f$y\f$.
     479             :  */
     480             : template <typename ArgumentTag>
     481           1 : struct ApplySwshJacobianInplace<Spectral::Swsh::Tags::Derivative<
     482             :     ArgumentTag, Spectral::Swsh::Tags::EthEth>> {
     483           0 :   using pre_swsh_derivative_tags = tmpl::list<>;
     484           0 :   using swsh_derivative_tags = tmpl::list<>;
     485           0 :   using integration_independent_tags =
     486             :       tmpl::list<Tags::OneMinusY, Tags::EthRDividedByR,
     487             :                  Tags::EthEthRDividedByR>;
     488             : 
     489           0 :   using return_tags = tmpl::list<Spectral::Swsh::Tags::Derivative<
     490             :       ArgumentTag, Spectral::Swsh::Tags::EthEth>>;
     491           0 :   using argument_tags = tmpl::append<integration_independent_tags>;
     492           0 :   using on_demand_argument_tags =
     493             :       tmpl::list<Tags::Dy<ArgumentTag>, Tags::Dy<Tags::Dy<ArgumentTag>>,
     494             :                  Spectral::Swsh::Tags::Derivative<Tags::Dy<ArgumentTag>,
     495             :                                                   Spectral::Swsh::Tags::Eth>>;
     496             : 
     497           0 :   static constexpr int spin =
     498             :       Spectral::Swsh::Tags::Derivative<ArgumentTag,
     499             :                                        Spectral::Swsh::Tags::EthEth>::spin;
     500             :   template <typename DyArgumentType, typename DyDyArgumentType,
     501             :             typename EthDyArgumentType>
     502           0 :   static void apply(
     503             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, spin>>*>
     504             :           eth_eth_argument,
     505             :       const Scalar<SpinWeighted<ComplexDataVector, 0>>& one_minus_y,
     506             :       const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_r_divided_by_r,
     507             :       const Scalar<SpinWeighted<ComplexDataVector, 2>>& eth_eth_r_divided_by_r,
     508             :       const DyArgumentType& dy_argument, const DyDyArgumentType& dy_dy_argument,
     509             :       const EthDyArgumentType& eth_dy_argument) {
     510             :     get(*eth_eth_argument) -=
     511             :         square(get(eth_r_divided_by_r)) *
     512             :             (square(get(one_minus_y)) *
     513             :              SpinWeighted<DyDyArgumentType, spin - 2>{dy_dy_argument}) +
     514             :         get(one_minus_y) *
     515             :             (2.0 * get(eth_r_divided_by_r) *
     516             :                  SpinWeighted<EthDyArgumentType, spin - 1>{eth_dy_argument} +
     517             :              get(eth_eth_r_divided_by_r) *
     518             :                  SpinWeighted<DyArgumentType, spin - 2>{dy_argument});
     519             :   }
     520             : };
     521             : 
     522             : /*!
     523             :  * \brief Specialization for the spin-weighted derivative \f$\bar{\eth}
     524             :  * \bar{\eth}\f$.
     525             :  *
     526             :  * \details The implemented equation is:
     527             :  *
     528             :  * \f[
     529             :  * \bar{\eth} \bar{\eth} F = \bar{\eth}^\prime \bar{\eth}^\prime F
     530             :  * - (1 - y)^2 \frac{(\bar{\eth} R)^2}{R^2} \partial_y^2 F
     531             :  * - (1 - y) \left( 2 \frac{\bar{\eth} R}{R} \bar{\eth} \partial_y F
     532             :  * + \frac{\bar{\eth} \bar{\eth} R}{R} \partial_y F\right),
     533             :  * \f]
     534             :  *
     535             :  * where \f$\bar{\eth} \bar{\eth}\f$ is the derivative at constant Bondi radius
     536             :  * \f$r\f$ and \f$\bar{\eth}^\prime \bar{\eth}^\prime\f$ is the derivative at
     537             :  * constant numerical radius \f$y\f$.
     538             :  */
     539             : template <typename ArgumentTag>
     540           1 : struct ApplySwshJacobianInplace<Spectral::Swsh::Tags::Derivative<
     541             :     ArgumentTag, Spectral::Swsh::Tags::EthbarEthbar>> {
     542           0 :   using pre_swsh_derivative_tags = tmpl::list<>;
     543           0 :   using swsh_derivative_tags = tmpl::list<>;
     544           0 :   using integration_independent_tags =
     545             :       tmpl::list<Tags::OneMinusY, Tags::EthRDividedByR,
     546             :                  Tags::EthEthRDividedByR>;
     547             : 
     548           0 :   using return_tags = tmpl::list<Spectral::Swsh::Tags::Derivative<
     549             :       ArgumentTag, Spectral::Swsh::Tags::EthbarEthbar>>;
     550           0 :   using argument_tags = tmpl::append<integration_independent_tags>;
     551           0 :   using on_demand_argument_tags =
     552             :       tmpl::list<Tags::Dy<ArgumentTag>, Tags::Dy<Tags::Dy<ArgumentTag>>,
     553             :                  Spectral::Swsh::Tags::Derivative<
     554             :                      Tags::Dy<ArgumentTag>, Spectral::Swsh::Tags::Ethbar>>;
     555             : 
     556           0 :   static constexpr int spin = Spectral::Swsh::Tags::Derivative<
     557             :       ArgumentTag, Spectral::Swsh::Tags::EthbarEthbar>::spin;
     558             :   template <typename DyArgumentType, typename DyDyArgumentType,
     559             :             typename EthbarDyArgumentType>
     560           0 :   static void apply(
     561             :       const gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, spin>>*>
     562             :           ethbar_ethbar_argument,
     563             :       const Scalar<SpinWeighted<ComplexDataVector, 0>>& one_minus_y,
     564             :       const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_r_divided_by_r,
     565             :       const Scalar<SpinWeighted<ComplexDataVector, 2>>& eth_eth_r_divided_by_r,
     566             :       const DyArgumentType& dy_argument, const DyDyArgumentType& dy_dy_argument,
     567             :       const EthbarDyArgumentType& ethbar_dy_argument) {
     568             :     get(*ethbar_ethbar_argument) -=
     569             :         square(conj(get(eth_r_divided_by_r))) *
     570             :             (square(get(one_minus_y)) *
     571             :              SpinWeighted<DyDyArgumentType, spin + 2>{dy_dy_argument}) +
     572             :         get(one_minus_y) *
     573             :             (2.0 * conj(get(eth_r_divided_by_r)) *
     574             :                  SpinWeighted<EthbarDyArgumentType, spin + 1>{
     575             :                      ethbar_dy_argument} +
     576             :              conj(get(eth_eth_r_divided_by_r)) *
     577             :                  SpinWeighted<DyArgumentType, spin + 2>{dy_argument});
     578             :   }
     579             : };
     580             : 
     581             : namespace detail {
     582             : // A helper to forward to the `ApplySwshJacobianInplace` mutators that takes
     583             : // advantage of the `OnDemandInputsForSwshJacobian`, which computes blaze
     584             : // template expressions for those quantities for which it is anticipated to be
     585             : // sufficiently cheap to repeatedly compute that it is worth saving the cost of
     586             : // additional storage (e.g. conjugates and derivatives for which the product
     587             : // rule applies)
     588             : template <typename DerivativeTag, typename DataBoxTagList,
     589             :           typename... OnDemandTags>
     590             : void apply_swsh_jacobian_helper(
     591             :     const gsl::not_null<db::DataBox<DataBoxTagList>*> box,
     592             :     tmpl::list<OnDemandTags...> /*meta*/) {
     593             :   db::mutate_apply<ApplySwshJacobianInplace<DerivativeTag>>(
     594             :       box, OnDemandInputsForSwshJacobian<OnDemandTags>{}(*box)...);
     595             : }
     596             : }  // namespace detail
     597             : 
     598             : /*!
     599             :  * \brief This routine evaluates the set of inputs to the CCE integrand for
     600             :  * `BondiValueTag` which are spin-weighted angular derivatives.
     601             : 
     602             :  * \details This function is called on the \ref DataBoxGroup holding the
     603             :  * relevant CCE data during each hypersurface integration step, after evaluating
     604             :  * `mutate_all_pre_swsh_derivatives_for_tag()` with template argument
     605             :  * `BondiValueTag` and before evaluating `ComputeBondiIntegrand<BondiValueTag>`.
     606             :  * Provided a \ref DataBoxGroup with the appropriate tags (including
     607             :  * `Cce::all_pre_swsh_derivative_tags`, `Cce::all_swsh_derivative_tags`,
     608             :  * `Cce::all_transform_buffer_tags`,  `Cce::pre_computation_tags`, and
     609             :  * `Cce::Tags::LMax`), this function will apply all of the necessary
     610             :  * mutations to update
     611             :  * `Cce::single_swsh_derivative_tags_to_compute_for<BondiValueTag>` and
     612             :  * `Cce::second_swsh_derivative_tags_to_compute_for<BondiValueTag>` to their
     613             :  * correct values for the current values of the remaining (input) tags.
     614             :  */
     615             : template <typename BondiValueTag, typename DataBoxTagList>
     616           1 : void mutate_all_swsh_derivatives_for_tag(
     617             :     const gsl::not_null<db::DataBox<DataBoxTagList>*> box) {
     618             :   // The collection of spin-weighted derivatives cannot be applied as individual
     619             :   // compute items, because it is better to aggregate similar spins and dispatch
     620             :   // to libsharp in groups. So, we supply a bulk mutate operation which takes in
     621             :   // multiple Variables from the presumed DataBox, and alters their values as
     622             :   // necessary.
     623             :   db::mutate_apply<Spectral::Swsh::AngularDerivatives<
     624             :       single_swsh_derivative_tags_to_compute_for_t<BondiValueTag>>>(box);
     625             :   tmpl::for_each<single_swsh_derivative_tags_to_compute_for_t<BondiValueTag>>(
     626             :       [&box](auto derivative_tag_v) {
     627             :         using derivative_tag = typename decltype(derivative_tag_v)::type;
     628             :         detail::apply_swsh_jacobian_helper<derivative_tag>(
     629             :             box, typename ApplySwshJacobianInplace<
     630             :                      derivative_tag>::on_demand_argument_tags{});
     631             :       });
     632             : 
     633             :   db::mutate_apply<Spectral::Swsh::AngularDerivatives<
     634             :       second_swsh_derivative_tags_to_compute_for_t<BondiValueTag>>>(box);
     635             :   tmpl::for_each<second_swsh_derivative_tags_to_compute_for_t<BondiValueTag>>(
     636             :       [&box](auto derivative_tag_v) {
     637             :         using derivative_tag = typename decltype(derivative_tag_v)::type;
     638             :         detail::apply_swsh_jacobian_helper<derivative_tag>(
     639             :             box, typename ApplySwshJacobianInplace<
     640             :                      derivative_tag>::on_demand_argument_tags{});
     641             :       });
     642             : }
     643             : }  // namespace Cce

Generated by: LCOV version 1.14