diff --git a/docs/References.bib b/docs/References.bib index 618ac4320ef2..2d980cc373cc 100644 --- a/docs/References.bib +++ b/docs/References.bib @@ -85,6 +85,21 @@ @article{Beckwith2011iy volume = "193", } +@article{Bishop1997ik, + author = "Bishop, Nigel T. and Gomez, Roberto and Lehner, Luis and + Maharaj, Manoj and Winicour, Jeffrey", + title = "{High powered gravitational news}", + journal = "Phys. Rev.", + volume = "D56", + year = "1997", + pages = "6298-6309", + doi = "10.1103/PhysRevD.56.6298", + eprint = "gr-qc/9708065", + archivePrefix = "arXiv", + primaryClass = "gr-qc", + url = "https://doi.org/10.1103/PhysRevD.56.6298", +} + @book{Cockburn1999, author = "Cockburn, Bernardo", editor = "Barth, Timothy J. and Deconinck, Herman", @@ -298,6 +313,22 @@ @article{Hemberger2012jz primaryClass = "gr-qc", } +@article{Handmer2014qha, + author = "Handmer, Casey J. and Szilagyi, B.", + title = "{Spectral Characteristic Evolution: A New Algorithm for + Gravitational Wave Propagation}", + journal = "Class. Quant. Grav.", + volume = "32", + year = "2015", + number = "2", + pages = "025008", + doi = "10.1088/0264-9381/32/2/025008", + eprint = "1406.7029", + archivePrefix = "arXiv", + primaryClass = "gr-qc", + url = "https://doi.org/10.1088/0264-9381/32/2/025008" +} + @book{Hartle2003gravity, author = "Hartle, J.B.", title = "Gravity: An Introduction to {Einstein's} General Relativity", diff --git a/src/DataStructures/Tags.hpp b/src/DataStructures/Tags.hpp index 280227114de4..2248b8d4231c 100644 --- a/src/DataStructures/Tags.hpp +++ b/src/DataStructures/Tags.hpp @@ -6,7 +6,9 @@ #include #include "DataStructures/DataBox/DataBoxTag.hpp" +#include "DataStructures/SpinWeighted.hpp" #include "DataStructures/Tensor/Metafunctions.hpp" +#include "DataStructures/Tensor/TypeAliases.hpp" /// \cond class DataVector; @@ -40,4 +42,52 @@ struct Modal : db::SimpleTag, db::PrefixTag { return "Modal(" + NodalTag::name() + ")"; } }; + +/// Given a Tag with a `type` of `Tensor`, acts as a new +/// version of the tag with `type` of `Tensor, ...>`, which is the preferred tensor type associated +/// with spin-weighted quantities. Here, `SpinConstant` must be a +/// `std::integral_constant` or similar type wrapper for a compile-time +/// constant, in order to work properly with DataBox utilities. +/// \note There are restrictions of which vector types can be wrapped +/// by SpinWeighted in a Tensor, so some `Tag`s that have valid `type`s may +/// give rise to compilation errors when used as `Tags::SpinWeighted`. If you find such trouble, consult the whitelist of possible +/// Tensor storage types in `Tensor.hpp`. +template +struct SpinWeighted : db::PrefixTag, db::SimpleTag { + static_assert(not is_any_spin_weighted_v, + "The SpinWeighted tag should only be used to create a " + "spin-weighted version of a non-spin-weighted tag. The " + "provided tag already has a spin-weighted type"); + using type = TensorMetafunctions::swap_type< + ::SpinWeighted, + typename Tag::type>; + using tag = Tag; + static std::string name() noexcept { + return "SpinWeighted(" + Tag::name() + ", " + + std::to_string(SpinConstant::value) + ")"; + } +}; + +namespace detail { +template +using product_t = + Scalar<::SpinWeighted>; +} // namespace detail + +/// A prefix tag representing the product of two other tags. Note that if +/// non-spin-weighted types are needed in this tag, the type alias +/// `detail::product_t` should be generalized to give a reasonable type for +/// those cases, or template specializations should be constructed for this tag. +template +struct Multiplies : db::PrefixTag, db::SimpleTag { + using type = detail::product_t, db::item_type>; + using tag = LhsTag; + static std::string name() noexcept { + return "Multiplies(" + LhsTag::name() + ", " + RhsTag::name() + ")"; + } +}; + } // namespace Tags diff --git a/src/DataStructures/Variables.hpp b/src/DataStructures/Variables.hpp index c4208794efa7..2f9027f0ed7a 100644 --- a/src/DataStructures/Variables.hpp +++ b/src/DataStructures/Variables.hpp @@ -909,7 +909,7 @@ struct Subitems -struct TempTensor { +struct TempTensor : db::SimpleTag { using type = T; static std::string name() noexcept { return std::string("TempTensor") + std::to_string(N); diff --git a/src/Evolution/Systems/CMakeLists.txt b/src/Evolution/Systems/CMakeLists.txt index 1138f0a84ac4..720ef03012ed 100644 --- a/src/Evolution/Systems/CMakeLists.txt +++ b/src/Evolution/Systems/CMakeLists.txt @@ -2,6 +2,7 @@ # See LICENSE.txt for details. add_subdirectory(Burgers) +add_subdirectory(Cce) add_subdirectory(CurvedScalarWave) add_subdirectory(GeneralizedHarmonic) add_subdirectory(GrMhd) diff --git a/src/Evolution/Systems/Cce/CMakeLists.txt b/src/Evolution/Systems/Cce/CMakeLists.txt new file mode 100644 index 000000000000..75ecfe6eb422 --- /dev/null +++ b/src/Evolution/Systems/Cce/CMakeLists.txt @@ -0,0 +1,16 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +set(LIBRARY Cce) + +set(LIBRARY_SOURCES + Equations.cpp + ) + +add_spectre_library(${LIBRARY} ${LIBRARY_SOURCES}) + +target_link_libraries( + ${LIBRARY} + INTERFACE DataStructures + INTERFACE ErrorHandling + ) diff --git a/src/Evolution/Systems/Cce/Equations.cpp b/src/Evolution/Systems/Cce/Equations.cpp new file mode 100644 index 000000000000..cad24f485915 --- /dev/null +++ b/src/Evolution/Systems/Cce/Equations.cpp @@ -0,0 +1,268 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "Evolution/Systems/Cce/Equations.hpp" + +#include + +#include "DataStructures/ComplexDataVector.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/SpinWeighted.hpp" +#include "Utilities/ConstantExpressions.hpp" + +// IWYU pragma: no_forward_declare Cce::Tags::BondiBeta +// IWYU pragma: no_forward_declare Cce::Tags::BondiH +// IWYU pragma: no_forward_declare Cce::Tags::BondiQ +// IWYU pragma: no_forward_declare Cce::Tags::BondiU +// IWYU pragma: no_forward_declare Cce::Tags::BondiW +// IWYU pragma: no_forward_declare Cce::Tags::Integrand +// IWYU pragma: no_forward_declare Cce::Tags::LinearFactor +// IWYU pragma: no_forward_declare Cce::Tags::LinearFactorForConjugate +// IWYU pragma: no_forward_declare Cce::Tags::PoleOfIntegrand +// IWYU pragma: no_forward_declare Cce::Tags::RegularIntegrand +// IWYU pragma: no_forward_declare SpinWeighted + +namespace Cce { +// suppresses doxygen problems with these functions +/// \cond + +void ComputeBondiIntegrand>::apply_impl( + const gsl::not_null*> integrand_for_beta, + const SpinWeighted& dy_j, + const SpinWeighted& j, + const SpinWeighted& one_minus_y) noexcept { + *integrand_for_beta = + 0.125 * one_minus_y * + (dy_j * conj(dy_j) - + 0.25 * square(j * conj(dy_j) + conj(j) * dy_j) / (1.0 + j * conj(j))); +} + +void ComputeBondiIntegrand>::apply_impl( + const gsl::not_null*> + pole_of_integrand_for_q, + const SpinWeighted& eth_beta) noexcept { + *pole_of_integrand_for_q = -4.0 * eth_beta; +} + +void ComputeBondiIntegrand>::apply_impl( + const gsl::not_null*> + regular_integrand_for_q, + const gsl::not_null*> script_aq, + const SpinWeighted& dy_beta, + const SpinWeighted& dy_j, + const SpinWeighted& j, + const SpinWeighted& eth_dy_beta, + const SpinWeighted& eth_j_jbar, + const SpinWeighted& eth_jbar_dy_j, + const SpinWeighted& ethbar_dy_j, + const SpinWeighted& ethbar_j, + const SpinWeighted& eth_r_divided_by_r, + const SpinWeighted& k) noexcept { + *script_aq = + 0.25 * (j * conj(ethbar_dy_j) - eth_jbar_dy_j - conj(ethbar_j) * dy_j + + 0.5 * eth_j_jbar * (conj(j) * dy_j + j * conj(dy_j)) / + (1.0 + j * conj(j)) - + (conj(j) * dy_j - j * conj(dy_j)) * eth_r_divided_by_r); + + *regular_integrand_for_q = + -2.0 * (*script_aq + j * conj(*script_aq) / k - eth_dy_beta + + 0.5 * ethbar_dy_j / k - dy_beta * eth_r_divided_by_r + + 0.5 * dy_j * conj(eth_r_divided_by_r) / k); +} + +void ComputeBondiIntegrand>::apply_impl( + const gsl::not_null*> + regular_integrand_for_u, + const SpinWeighted& exp_2_beta, + const SpinWeighted& j, + const SpinWeighted& q, + const SpinWeighted& k, + const SpinWeighted& r) noexcept { + *regular_integrand_for_u = 0.5 * exp_2_beta / r * (k * q - j * conj(q)); +} + +void ComputeBondiIntegrand>::apply_impl( + const gsl::not_null*> + pole_of_integrand_for_w, + const SpinWeighted& ethbar_u) noexcept { + *pole_of_integrand_for_w = ethbar_u + conj(ethbar_u); +} + +void ComputeBondiIntegrand>::apply_impl( + const gsl::not_null*> + regular_integrand_for_w, + const gsl::not_null*> script_av, + const SpinWeighted& dy_u, + const SpinWeighted& exp_2_beta, + const SpinWeighted& j, + const SpinWeighted& q, + const SpinWeighted& eth_beta, + const SpinWeighted& eth_eth_beta, + const SpinWeighted& eth_ethbar_beta, + const SpinWeighted& eth_ethbar_j, + const SpinWeighted& eth_ethbar_j_jbar, + const SpinWeighted& eth_j_jbar, + const SpinWeighted& ethbar_dy_u, + const SpinWeighted& ethbar_ethbar_j, + const SpinWeighted& ethbar_j, + const SpinWeighted& eth_r_divided_by_r, + const SpinWeighted& k, + const SpinWeighted& r) noexcept { + *script_av = + (eth_beta * conj(ethbar_j) + 0.5 * ethbar_ethbar_j + + j * square(conj(eth_beta)) + j * conj(eth_eth_beta) + + 0.125 * eth_j_jbar * conj(eth_j_jbar) / (k * (1.0 + j * conj(j))) + + 0.5 * + (1.0 - 0.25 * eth_ethbar_j_jbar - eth_j_jbar * conj(eth_beta) - + 0.5 * conj(ethbar_j) * ethbar_j - 0.5 * conj(j) * eth_ethbar_j) / + k + + k * (0.5 - eth_ethbar_beta - eth_beta * conj(eth_beta) - + 0.25 * q * conj(q)) + + 0.25 * j * square(conj(q))); + + *regular_integrand_for_w = + 0.5 * (0.5 * (ethbar_dy_u + conj(ethbar_dy_u) + + conj(dy_u) * eth_r_divided_by_r + + dy_u * conj(eth_r_divided_by_r)) - + 1.0 / r + 0.5 * exp_2_beta * (*script_av + conj(*script_av)) / r); +} + +void ComputeBondiIntegrand>::apply_impl( + const gsl::not_null*> + pole_of_integrand_for_h, + const SpinWeighted& j, + const SpinWeighted& u, + const SpinWeighted& w, + const SpinWeighted& eth_u, + const SpinWeighted& ethbar_j, + const SpinWeighted& ethbar_jbar_u, + const SpinWeighted& ethbar_u, + const SpinWeighted& k) noexcept { + *pole_of_integrand_for_h = -0.5 * conj(ethbar_jbar_u) - j * conj(ethbar_u) - + 0.5 * j * ethbar_u - k * eth_u - + 0.5 * u * ethbar_j + 2.0 * j * w; +} + +void ComputeBondiIntegrand>::apply_impl( + const gsl::not_null*> + regular_integrand_for_h, + const gsl::not_null*> script_aj, + const gsl::not_null*> script_bj, + const gsl::not_null*> script_cj, + const SpinWeighted& dy_dy_j, + const SpinWeighted& dy_j, + const SpinWeighted& dy_w, + const SpinWeighted& exp_2_beta, + const SpinWeighted& j, + const SpinWeighted& q, + const SpinWeighted& u, + const SpinWeighted& w, + const SpinWeighted& eth_beta, + const SpinWeighted& eth_eth_beta, + const SpinWeighted& eth_ethbar_beta, + const SpinWeighted& eth_ethbar_j, + const SpinWeighted& eth_ethbar_j_jbar, + const SpinWeighted& eth_j_jbar, + const SpinWeighted& eth_q, + const SpinWeighted& eth_u, + const SpinWeighted& eth_ubar_dy_j, + const SpinWeighted& ethbar_dy_j, + const SpinWeighted& ethbar_ethbar_j, + const SpinWeighted& ethbar_j, + const SpinWeighted& ethbar_jbar_dy_j, + const SpinWeighted& ethbar_jbar_q_minus_2_eth_beta, + const SpinWeighted& ethbar_q, + const SpinWeighted& ethbar_u, + const SpinWeighted& du_r_divided_by_r, + const SpinWeighted& eth_r_divided_by_r, + const SpinWeighted& k, + const SpinWeighted& one_minus_y, + const SpinWeighted& r) noexcept { + *script_aj = + 0.25 * + (conj(ethbar_ethbar_j) - + 0.25 * (4.0 + eth_ethbar_j_jbar - j * conj(eth_ethbar_j)) / + (k * (1.0 + j * conj(j))) + + (3.0 - eth_ethbar_beta - + conj(j) * eth_ethbar_j * (1.0 - 0.25 / (1.0 + j * conj(j)))) / + k + + conj(ethbar_j) * (2.0 * eth_beta + + 0.5 * + (j * conj(eth_j_jbar) - + ethbar_j * (2.0 * (1.0 + j * conj(j)) - 1.0)) / + (k * (1.0 + j * conj(j))) - + q)); + + *script_bj = + 0.25 * (2.0 * dy_w - + conj(j) * eth_u * (conj(j) * dy_j + j * conj(dy_j)) / k + + 1.0 / r + u * ethbar_j * conj(dy_j) - + 0.5 * u * conj(eth_j_jbar) * (conj(j) * dy_j + j * conj(dy_j)) / + (1.0 + j * conj(j)) + + conj(u) * (conj(ethbar_jbar_dy_j) - j * conj(ethbar_dy_j))) + + one_minus_y * 0.25 * + (du_r_divided_by_r * dy_j * + (conj(j) * (conj(j) * dy_j + j * conj(dy_j)) / + (1.0 + j * conj(j)) - + 2.0 * conj(dy_j)) - + w * (dy_j * conj(dy_j) - + 0.25 * square((conj(j) * dy_j + j * conj(dy_j))) / + (1.0 + j * conj(j)))) + + square(one_minus_y) * 0.125 * + (0.25 * square(conj(j) * dy_j + j * conj(dy_j)) / + (1.0 + j * conj(j)) - + dy_j * conj(dy_j)) / + r; + + *script_cj = 0.5 * ethbar_j * k * (eth_beta - 0.5 * q); + + *regular_integrand_for_h = + j * (*script_bj + conj(*script_bj)) - + 0.5 * (eth_ubar_dy_j + u * ethbar_dy_j + + u * dy_j * conj(eth_r_divided_by_r)) + + 0.5 * exp_2_beta / r * + (*script_cj + square(j) / (1.0 + j * conj(j)) * conj(*script_cj) - + j * (*script_aj + conj(*script_aj)) + eth_eth_beta - 0.5 * eth_q + + 0.25 * (conj(ethbar_jbar_q_minus_2_eth_beta) - j * conj(ethbar_q)) / + k + + square(eth_beta - 0.5 * q)) - + dy_j * 0.5 * + (conj(j) * eth_u / k - j * k * conj(eth_u) + + 0.5 * (1.0 + j * conj(j)) * (conj(ethbar_u) - ethbar_u) + + conj(u) * eth_r_divided_by_r - w) + + conj(dy_j) * (0.5 * j * eth_u * (j * conj(j) / k) - + 0.25 * square(j) * (ethbar_u - conj(ethbar_u))) + + one_minus_y * + (0.5 * (dy_dy_j * (w + 2.0 * du_r_divided_by_r) - dy_j / r) + + 0.5 * dy_j * (dy_w + 1.0 / r)) + + square(one_minus_y) * 0.25 * dy_dy_j / r; +} + +void ComputeBondiIntegrand>::apply_impl( + const gsl::not_null*> + linear_factor_for_h, + const gsl::not_null*> script_djbar, + const SpinWeighted& dy_j, + const SpinWeighted& j, + const SpinWeighted& one_minus_y) noexcept { + *script_djbar = 0.25 * one_minus_y * + (-2.0 * dy_j + + j * (conj(j) * dy_j + j * conj(dy_j)) / (1.0 + j * conj(j))); + *linear_factor_for_h = 1.0 + j * conj(*script_djbar); +} + +void ComputeBondiIntegrand>:: + apply_impl( + const gsl::not_null*> + linear_factor_for_conjugate_h, + const gsl::not_null*> script_djbar, + const SpinWeighted& dy_j, + const SpinWeighted& j, + const SpinWeighted& one_minus_y) noexcept { + *script_djbar = 0.25 * one_minus_y * + (-2.0 * dy_j + + j * (conj(j) * dy_j + j * conj(dy_j)) / (1.0 + j * conj(j))); + *linear_factor_for_conjugate_h = j * (*script_djbar); +} +/// \endcond +} // namespace Cce diff --git a/src/Evolution/Systems/Cce/Equations.hpp b/src/Evolution/Systems/Cce/Equations.hpp new file mode 100644 index 000000000000..f86f317fab63 --- /dev/null +++ b/src/Evolution/Systems/Cce/Equations.hpp @@ -0,0 +1,959 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include + +#include "DataStructures/SpinWeighted.hpp" // IWYU pragma: keep +#include "DataStructures/Tags.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Variables.hpp" +#include "Evolution/Systems/Cce/Tags.hpp" // IWYU pragma: keep +#include "NumericalAlgorithms/Spectral/SwshTags.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/TMPL.hpp" + +// IWYU pragma: no_forward_declare Cce::Tags::BondiBeta +// IWYU pragma: no_forward_declare Cce::Tags::DuRDividedByR +// IWYU pragma: no_forward_declare Cce::Tags::EthRDividedByR +// IWYU pragma: no_forward_declare Cce::Tags::Exp2Beta +// IWYU pragma: no_forward_declare Cce::Tags::BondiH +// IWYU pragma: no_forward_declare Cce::Tags::BondiJ +// IWYU pragma: no_forward_declare Cce::Tags::BondiJbar +// IWYU pragma: no_forward_declare Cce::Tags::JbarQMinus2EthBeta +// IWYU pragma: no_forward_declare Cce::Tags::BondiK +// IWYU pragma: no_forward_declare Cce::Tags::OneMinusY +// IWYU pragma: no_forward_declare Cce::Tags::BondiQ +// IWYU pragma: no_forward_declare Cce::Tags::BondiR +// IWYU pragma: no_forward_declare Cce::Tags::BondiU +// IWYU pragma: no_forward_declare Cce::Tags::BondiUbar +// IWYU pragma: no_forward_declare Cce::Tags::BondiW +// IWYU pragma: no_forward_declare ::Tags::Multiplies +// IWYU pragma: no_forward_declare Cce::Tags::Dy +// IWYU pragma: no_forward_declare Cce::Tags::Integrand +// IWYU pragma: no_forward_declare Cce::Tags::LinearFactor +// IWYU pragma: no_forward_declare Cce::Tags::LinearFactorForConjugate +// IWYU pragma: no_forward_declare Cce::Tags::PoleOfIntegrand +// IWYU pragma: no_forward_declare Cce::Tags::RegularIntegrand +// IWYU pragma: no_forward_declare Spectral::Swsh::Tags::Eth +// IWYU pragma: no_forward_declare Spectral::Swsh::Tags::EthEth +// IWYU pragma: no_forward_declare Spectral::Swsh::Tags::EthEthbar +// IWYU pragma: no_forward_declare Spectral::Swsh::Tags::Ethbar +// IWYU pragma: no_forward_declare Spectral::Swsh::Tags::EthbarEthbar +// IWYU pragma: no_forward_declare Spectral::Swsh::Tags::Derivative +// IWYU pragma: no_forward_declare Tags::TempTensor +// IWYU pragma: no_forward_declare Tags::SpinWeighted +// IWYU pragma: no_forward_declare SpinWeighted +// IWYU pragma: no_forward_declare Tensor + +/// \cond +class ComplexDataVector; +/// \endcond + +/// Contains functionality for Cauchy Characteristic Extraction +namespace Cce { + +namespace detail { +template +struct integrand_terms_to_compute_for_bondi_variable_impl; + +// template specializations for the individual tags +template <> +struct integrand_terms_to_compute_for_bondi_variable_impl { + using type = tmpl::list>; +}; +template <> +struct integrand_terms_to_compute_for_bondi_variable_impl { + using type = tmpl::list, + Tags::RegularIntegrand>; +}; +template <> +struct integrand_terms_to_compute_for_bondi_variable_impl { + using type = tmpl::list>; +}; +template <> +struct integrand_terms_to_compute_for_bondi_variable_impl { + using type = tmpl::list, + Tags::RegularIntegrand>; +}; +template <> +struct integrand_terms_to_compute_for_bondi_variable_impl { + using type = tmpl::list, + Tags::RegularIntegrand, + Tags::LinearFactor, + Tags::LinearFactorForConjugate>; +}; +} // namespace detail + +/// \brief A struct for providing a `tmpl::list` of integrand tags that need to +/// be computed before integration can proceed for a given Bondi variable tag. +template +using integrand_terms_to_compute_for_bondi_variable = + typename detail::integrand_terms_to_compute_for_bondi_variable_impl< + BondiVariable>::type; + +/*! + * \brief Computes one of the inputs for the integration of one of the + * Characteristic hypersurface equations. + * + * \details The template argument must be one of the integrand-related prefix + * tags templated on a Bondi quantity tag for which that integrand is required. + * The relevant prefix tags are `Tags::Integrand`, `Tags::PoleOfIntegrand`, + * `Tags::RegularIntegrand`, `Tags::LinearFactor`, and + * `Tags::LinearFactorForConjugate`. The Bondi quantity tags that these tags may + * wrap are `Tags::BondiBeta`, `Tags::BondiQ`, `Tags::BondiU`, `Tags::BondiW`, + * and `Tags::BondiH`. + * + * The integrand terms which may be computed for a given Bondi variable are + * enumerated in the type alias `integrand_terms_to_compute_for_bondi_variable`, + * which takes as a single template argument the tag for which integrand terms + * would be computed, and is a `tmpl::list` of the integrand terms needed. + * + * The resulting quantity is returned by `not_null` pointer, and the required + * argument tags are given in `return_tags` and `argument_tags` type aliases, + * where the `return_tags` are passed by `not_null` pointer (so include + * temporary buffers) and the `argument_tags` are passed by const reference. + * + * Additional mathematical details for each of the computations can be found in + * the template specializations of this struct. All of the specializations have + * a static `apply` function which takes as arguments + * `SpinWeighted`s for each of the Bondi arguments. + */ +template +struct ComputeBondiIntegrand; + +/*! + * \brief Computes the integrand (right-hand side) of the equation which + * determines the radial (y) dependence of the Bondi quantity \f$\beta\f$. + * + * \details The quantity \f$\beta\f$ is defined via the Bondi form of the + * metric: + * \f[ + * ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2 + * 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] + * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical + * angular dyad \f$q^A\f$: + * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A + * \bar{q}^B.\f] + * See \cite Bishop1997ik \cite Handmer2014qha for full details. + * + * We write the equations of motion in the compactified coordinate \f$ y \equiv + * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the + * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of + * the worldtube. The equation which determines \f$\beta\f$ on a surface of + * constant \f$u\f$ given \f$J\f$ on the same surface is + * \f[\partial_y (\beta) = + * \frac{1}{8} (-1 + y) \left(\partial_y (J) \partial_y(\bar{J}) + * - \frac{(\partial_y (J \bar{J}))^2}{4 K^2}\right). \f] + */ +template <> +struct ComputeBondiIntegrand> { + public: + using pre_swsh_derivative_tags = + tmpl::list, Tags::BondiJ>; + using swsh_derivative_tags = tmpl::list<>; + using integration_independent_tags = tmpl::list; + using temporary_tags = tmpl::list<>; + + using return_tags = tmpl::append>, + temporary_tags>; + using argument_tags = + tmpl::append; + + template + static void apply( + const gsl::not_null>*> + integrand_for_beta, + const Args&... args) noexcept { + apply_impl(make_not_null(&get(*integrand_for_beta)), get(args)...); + } + + private: + static void apply_impl( + gsl::not_null*> integrand_for_beta, + const SpinWeighted& dy_j, + const SpinWeighted& j, + const SpinWeighted& one_minus_y) noexcept; +}; + +/*! + * \brief Computes the pole part of the integrand (right-hand side) of the + * equation which determines the radial (y) dependence of the Bondi quantity + * \f$Q\f$. + * + * \details The quantity \f$Q\f$ is defined via the Bondi form of the metric: + * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2 + * 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] + * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical + * angular dyad \f$q^A\f$: + * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B,\f] + * and \f$Q\f$ is defined as a supplemental variable for radial integration of + * \f$U\f$: + * \f[ Q_A = r^2 e^{-2\beta} h_{AB} \partial_r U^B\f] + * and \f$Q = Q_A q^A\f$. See \cite Bishop1997ik \cite Handmer2014qha for full + * details. + * + * We write the equations of motion in the compactified coordinate \f$ y \equiv + * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the + * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of + * the worldtube. The equation which determines \f$Q\f$ on a surface of constant + * \f$u\f$ given \f$J\f$ and \f$\beta\f$ on the same surface is written as + * \f[(1 - y) \partial_y Q + 2 Q = A_Q + (1 - y) B_Q.\f] + * We refer to \f$A_Q\f$ as the "pole part" of the integrand and \f$B_Q\f$ + * as the "regular part". The pole part is computed by this function, and has + * the expression + * \f[A_Q = -4 \eth \beta.\f] + */ +template <> +struct ComputeBondiIntegrand> { + public: + using pre_swsh_derivative_tags = tmpl::list<>; + using swsh_derivative_tags = + tmpl::list>; + using integration_independent_tags = tmpl::list<>; + using temporary_tags = tmpl::list<>; + + using return_tags = + tmpl::append>, + temporary_tags>; + using argument_tags = + tmpl::append; + + template + static void apply( + const gsl::not_null>*> + pole_of_integrand_for_q, + const Args&... args) noexcept { + apply_impl(make_not_null(&get(*pole_of_integrand_for_q)), get(args)...); + } + + private: + static void apply_impl( + gsl::not_null*> + pole_of_integrand_for_q, + const SpinWeighted& eth_beta) noexcept; +}; + +/*! + * \brief Computes the regular part of the integrand (right-hand side) of the + * equation which determines the radial (y) dependence of the Bondi quantity + * \f$Q\f$. + * + * \details The quantity \f$Q\f$ is defined via the Bondi form of the metric: + * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2 + * 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] + * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical + * angular dyad \f$q^A\f$: + * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B,\f] + * and \f$Q\f$ is defined as a supplemental variable for radial integration of + * \f$U\f$: + * \f[ Q_A = r^2 e^{-2\beta} h_{AB} \partial_r U^B\f] + * and \f$Q = Q_A q^A\f$. See \cite Bishop1997ik \cite Handmer2014qha for + * full details. + * + * We write the equations of motion in the compactified coordinate \f$ y \equiv + * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the + * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of + * the worldtube. The equation which determines \f$Q\f$ on a surface of constant + * \f$u\f$ given \f$J\f$ and \f$\beta\f$ on the same surface is written as + * \f[(1 - y) \partial_y Q + 2 Q = A_Q + (1 - y) B_Q. \f] + * We refer to \f$A_Q\f$ as the "pole part" of the integrand and \f$B_Q\f$ as + * the "regular part". The regular part is computed by this function, and has + * the expression + * \f[ B_Q = - \left(2 \mathcal{A}_Q + \frac{2 + * \bar{\mathcal{A}_Q} J}{K} - 2 \partial_y (\eth (\beta)) + + * \frac{\partial_y (\bar{\eth} (J))}{K}\right), \f] + * where + * \f[ \mathcal{A}_Q = - \tfrac{1}{4} \eth (\bar{J} \partial_y (J)) + + * \tfrac{1}{4} J \partial_y (\eth (\bar{J})) - \tfrac{1}{4} \eth (\bar{J}) + * \partial_y (J) + \frac{\eth (J \bar{J}) \partial_y (J \bar{J})}{8 K^2} - + * \frac{\bar{J} \eth (R) \partial_y (J)}{4 R}. \f]. + */ +template <> +struct ComputeBondiIntegrand> { + public: + using pre_swsh_derivative_tags = + tmpl::list, Tags::Dy, + Tags::BondiJ>; + using swsh_derivative_tags = tmpl::list< + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Eth>, + Spectral::Swsh::Tags::Derivative< + ::Tags::Multiplies, + Spectral::Swsh::Tags::Eth>, + Spectral::Swsh::Tags::Derivative< + ::Tags::Multiplies>, + Spectral::Swsh::Tags::Eth>, + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Ethbar>, + Spectral::Swsh::Tags::Derivative>; + using integration_independent_tags = + tmpl::list; + using temporary_tags = + tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>, + std::integral_constant>>; + + using return_tags = + tmpl::append>, + temporary_tags>; + using argument_tags = + tmpl::append; + + template + static void apply( + const gsl::not_null>*> + regular_integrand_for_q, + const gsl::not_null>*> + script_aq, + const Args&... args) noexcept { + apply_impl(make_not_null(&get(*regular_integrand_for_q)), + make_not_null(&get(*script_aq)), get(args)...); + } + + private: + static void apply_impl( + gsl::not_null*> + regular_integrand_for_q, + gsl::not_null*> script_aq, + const SpinWeighted& dy_beta, + const SpinWeighted& dy_j, + const SpinWeighted& j, + const SpinWeighted& eth_dy_beta, + const SpinWeighted& eth_j_jbar, + const SpinWeighted& eth_jbar_dy_j, + const SpinWeighted& ethbar_dy_j, + const SpinWeighted& ethbar_j, + const SpinWeighted& eth_r_divided_by_r, + const SpinWeighted& k) noexcept; +}; + +/*! + * \brief Computes the integrand (right-hand side) of the equation which + * determines the radial (y) dependence of the Bondi quantity \f$U\f$. + * + * \details The quantity \f$U\f$ is defined via the Bondi form of the metric: + * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2 + * 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] + * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical + * angular dyad \f$q^A\f$: + * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B,\f] + * and \f$Q\f$ is defined as a supplemental variable for radial integration of + * \f$U\f$: + * \f[ Q_A = r^2 e^{-2\beta} h_{AB} \partial_r U^B\f] + * and \f$U = U_A q^A\f$. See \cite Bishop1997ik \cite Handmer2014qha for full + * details. + * + * We write the equations of motion in the compactified coordinate \f$ y \equiv + * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the + * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of + * the worldtube. The equation which determines \f$U\f$ on a surface of constant + * \f$u\f$ given \f$J\f$, \f$\beta\f$, and \f$Q\f$ on the same surface is + * written as + * \f[\partial_y U = \frac{e^{2\beta}}{2 R} (K Q - J \bar{Q}). \f] + */ +template <> +struct ComputeBondiIntegrand> { + public: + using pre_swsh_derivative_tags = + tmpl::list; + using swsh_derivative_tags = tmpl::list<>; + using integration_independent_tags = tmpl::list; + using temporary_tags = tmpl::list<>; + + using return_tags = + tmpl::append>, temporary_tags>; + using argument_tags = + tmpl::append; + + template + static void apply( + const gsl::not_null>*> + regular_integrand_for_u, + const Args&... args) noexcept { + apply_impl(make_not_null(&get(*regular_integrand_for_u)), get(args)...); + } + + private: + static void apply_impl(gsl::not_null*> + regular_integrand_for_u, + const SpinWeighted& exp_2_beta, + const SpinWeighted& j, + const SpinWeighted& q, + const SpinWeighted& k, + const SpinWeighted& r) noexcept; +}; + +/*! + * \brief Computes the pole part of the integrand (right-hand side) of the + * equation which determines the radial (y) dependence of the Bondi quantity + * \f$W\f$. + * + * \details The quantity \f$W\f$ is defined via the Bondi form of the metric: + * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2 + * 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] + * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical + * angular dyad \f$q^A\f$: + * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B.\f] + * See \cite Bishop1997ik \cite Handmer2014qha for full details. + * + * We write the equations of motion in the compactified coordinate \f$ y \equiv + * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the + * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of + * the worldtube. The equation which determines \f$W\f$ on a surface of constant + * \f$u\f$ given \f$J\f$,\f$\beta\f$, \f$Q\f$, and \f$U\f$ on the same surface + * is written as + * \f[(1 - y) \partial_y W + 2 W = A_W + (1 - y) B_W.\f] We refer + * to \f$A_W\f$ as the "pole part" of the integrand and \f$B_W\f$ as the + * "regular part". The pole part is computed by this function, and has the + * expression + * \f[A_W = \eth (\bar{U}) + \bar{\eth} (U).\f] + */ +template <> +struct ComputeBondiIntegrand> { + public: + using pre_swsh_derivative_tags = tmpl::list<>; + using swsh_derivative_tags = tmpl::list>; + using integration_independent_tags = tmpl::list<>; + using temporary_tags = tmpl::list<>; + + using return_tags = + tmpl::append>, + temporary_tags>; + using argument_tags = + tmpl::append; + + template + static void apply( + const gsl::not_null>*> + pole_of_integrand_for_w, + const Args&... args) noexcept { + apply_impl(make_not_null(&get(*pole_of_integrand_for_w)), get(args)...); + } + + private: + static void apply_impl( + gsl::not_null*> + pole_of_integrand_for_w, + const SpinWeighted& ethbar_u) noexcept; +}; + +/*! + * \brief Computes the regular part of the integrand (right-hand side) of the + * equation which determines the radial (y) dependence of the Bondi quantity + * \f$W\f$. + * + * \details The quantity \f$W\f$ is defined via the Bondi form of the metric: + * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2 + * 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] + * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical + * angular dyad \f$q^A\f$: + * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B,\f] + * See \cite Bishop1997ik \cite Handmer2014qha for full details. + * + * We write the equations of motion in the compactified coordinate \f$ y \equiv + * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the + * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of + * the worldtube. The equation which determines \f$W\f$ on a surface of constant + * \f$u\f$ given \f$J\f$, \f$\beta\f$, \f$Q\f$, \f$U\f$ on the same surface is + * written as + * \f[(1 - y) \partial_y W + 2 W = A_W + (1 - y) B_W. \f] + * We refer to \f$A_W\f$ as the "pole part" of the integrand and \f$B_W\f$ as + * the "regular part". The regular part is computed by this function, and has + * the expression + * \f[ B_W = \tfrac{1}{4} \partial_y (\eth (\bar{U})) + \tfrac{1}{4} \partial_y + * (\bar{\eth} (U)) - \frac{1}{2 R} + \frac{e^{2 \beta} (\mathcal{A}_W + + * \bar{\mathcal{A}_W})}{4 R}, \f] + * where + * \f{align*} + * \mathcal{A}_W =& - \eth (\beta) \eth (\bar{J}) + \tfrac{1}{2} \bar{\eth} + * (\bar{\eth} (J)) + 2 \bar{\eth} (\beta) \bar{\eth} (J) + (\bar{\eth} + * (\beta))^2 J + \bar{\eth} + * (\bar{\eth} (\beta)) J + \frac{\eth (J \bar{J}) \bar{\eth} (J \bar{J})}{8 + * K^3} + \frac{1}{2 K} - \frac{\eth (\bar{\eth} (J \bar{J}))}{8 K} - + * \frac{\eth (J + * \bar{J}) \bar{\eth} (\beta)}{2 K} \nonumber \\ + * &- \frac{\eth (\bar{J}) \bar{\eth} (J)}{4 K} - \frac{\eth (\bar{\eth} (J)) + * \bar{J}}{4 K} + \tfrac{1}{2} K - \eth (\bar{\eth} (\beta)) K - \eth + * (\beta) \bar{\eth} (\beta) K + \tfrac{1}{4} (- K Q \bar{Q} + J \bar{Q}^2). + * \f} + */ +template <> +struct ComputeBondiIntegrand> { + public: + using pre_swsh_derivative_tags = + tmpl::list, Tags::Exp2Beta, Tags::BondiJ, + Tags::BondiQ>; + using swsh_derivative_tags = tmpl::list< + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Derivative< + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Eth>, + Spectral::Swsh::Tags::Derivative< + ::Tags::Multiplies, + Spectral::Swsh::Tags::EthEthbar>, + Spectral::Swsh::Tags::Derivative< + ::Tags::Multiplies, + Spectral::Swsh::Tags::Eth>, + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Ethbar>, + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Derivative>; + using integration_independent_tags = + tmpl::list; + using temporary_tags = + tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>, + std::integral_constant>>; + + using return_tags = + tmpl::append>, + temporary_tags>; + using argument_tags = + tmpl::append; + + template + static void apply( + const gsl::not_null>*> + regular_integrand_for_w, + const gsl::not_null>*> + script_av, + const Args&... args) noexcept { + apply_impl(make_not_null(&get(*regular_integrand_for_w)), + make_not_null(&get(*script_av)), get(args)...); + } + + private: + static void apply_impl( + gsl::not_null*> + regular_integrand_for_w, + gsl::not_null*> script_av, + const SpinWeighted& dy_u, + const SpinWeighted& exp_2_beta, + const SpinWeighted& j, + const SpinWeighted& q, + const SpinWeighted& eth_beta, + const SpinWeighted& eth_eth_beta, + const SpinWeighted& eth_ethbar_beta, + const SpinWeighted& eth_ethbar_j, + const SpinWeighted& eth_ethbar_j_jbar, + const SpinWeighted& eth_j_jbar, + const SpinWeighted& ethbar_dy_u, + const SpinWeighted& ethbar_ethbar_j, + const SpinWeighted& ethbar_j, + const SpinWeighted& eth_r_divided_by_r, + const SpinWeighted& k, + const SpinWeighted& r) noexcept; +}; + +/*! + * \brief Computes the pole part of the integrand (right-hand side) of the + * equation which determines the radial (y) dependence of the Bondi quantity + * \f$H\f$. + * + * \details The quantity \f$H \equiv \partial_u J\f$ (evaluated at constant y) + * is defined via the Bondi form of the metric: + * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2 + * 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] + * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical + * angular dyad \f$q^A\f$: + * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B.\f] + * See \cite Bishop1997ik \cite Handmer2014qha for full details. + * + * We write the equations of motion in the compactified coordinate \f$ y \equiv + * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the + * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of + * the worldtube. The equation which determines \f$W\f$ on a surface of constant + * \f$u\f$ given \f$J\f$,\f$\beta\f$, \f$Q\f$, \f$U\f$, and \f$W\f$ on the same + * surface is written as + * \f[(1 - y) \partial_y H + H + (1 - y)(\mathcal{D}_J H + * + \bar{\mathcal{D}}_J \bar{H}) = A_J + (1 - y) B_J.\f] + * + * We refer to \f$A_J\f$ as the "pole part" of the integrand + * and \f$B_J\f$ as the "regular part". The pole part is computed by this + * function, and has the expression + * \f{align*} + * A_J =& - \tfrac{1}{2} \eth (J \bar{U}) - \eth (\bar{U}) J - \tfrac{1}{2} + * \bar{\eth} (U) J - \eth (U) K - \tfrac{1}{2} (\bar{\eth} (J) U) + 2 J W + * \f} + */ +template <> +struct ComputeBondiIntegrand> { + public: + using pre_swsh_derivative_tags = + tmpl::list; + using swsh_derivative_tags = tmpl::list< + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Derivative< + ::Tags::Multiplies, + Spectral::Swsh::Tags::Ethbar>, + Spectral::Swsh::Tags::Derivative>; + using integration_independent_tags = tmpl::list; + using temporary_tags = tmpl::list<>; + + using return_tags = + tmpl::append>, + temporary_tags>; + using argument_tags = + tmpl::append; + + template + static void apply( + const gsl::not_null>*> + pole_of_integrand_for_h, + const Args&... args) noexcept { + apply_impl(make_not_null(&get(*pole_of_integrand_for_h)), get(args)...); + } + + private: + static void apply_impl( + gsl::not_null*> + pole_of_integrand_for_h, + const SpinWeighted& j, + const SpinWeighted& u, + const SpinWeighted& w, + const SpinWeighted& eth_u, + const SpinWeighted& ethbar_j, + const SpinWeighted& ethbar_jbar_u, + const SpinWeighted& ethbar_u, + const SpinWeighted& k) noexcept; +}; + +/*! + * \brief Computes the pole part of the integrand (right-hand side) of the + * equation which determines the radial (y) dependence of the Bondi quantity + * \f$H\f$. + * + * \details The quantity \f$H \equiv \partial_u J\f$ (evaluated at constant y) + * is defined via the Bondi form of the metric: + * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2 + * 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] + * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical + * angular dyad \f$q^A\f$: + * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B.\f] + * See \cite Bishop1997ik \cite Handmer2014qha for full details. + * + * We write the equations of motion in the compactified coordinate \f$ y \equiv + * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the + * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of + * the worldtube. The equation which determines \f$H\f$ on a surface of constant + * \f$u\f$ given \f$J\f$,\f$\beta\f$, \f$Q\f$, \f$U\f$, and \f$W\f$ on the same + * surface is written as + * \f[(1 - y) \partial_y H + H + (1 - y)(\mathcal{D}_J H + \bar{\mathcal{D}}_J + * \bar{H}) = A_J + (1 - y) B_J.\f] + * We refer to \f$A_J\f$ as the "pole part" of the integrand + * and \f$B_J\f$ as the "regular part". The pole part is computed by this + * function, and has the expression + * \f{align*} + * B_J =& -\tfrac{1}{2} \left(\eth(\partial_y (J) \bar{U}) + \partial_y + * (\bar{\eth} (J)) U \right) + J (\mathcal{B}_J + \bar{\mathcal{B}}_J) \notag\\ + * &+ \frac{e^{2 \beta}}{2 R} \left(\mathcal{C}_J + \eth (\eth (\beta)) - + * \tfrac{1}{2} \eth (Q) - (\mathcal{A}_J + \bar{\mathcal{A}_J}) J + + * \frac{\bar{\mathcal{C}_J} J^2}{K^2} + \frac{\eth (J (-2 \bar{\eth} (\beta) + + * \bar{Q}))}{4 K} - \frac{\eth (\bar{Q}) J}{4 K} + (\eth (\beta) - + * \tfrac{1}{2} Q)^2\right) \notag\\ + * &- \partial_y (J) \left(\frac{\eth (U) \bar{J}}{2 K} - \tfrac{1}{2} + * \bar{\eth} (\bar{U}) J K + \tfrac{1}{4} (\eth (\bar{U}) - \bar{\eth} (U)) + * K^2 + \frac{1}{2} \frac{\eth (R) \bar{U}}{R} - \frac{1}{2} + * W\right)\notag\\ + * &+ \partial_y (\bar{J}) \left(- \tfrac{1}{4} (- \eth (\bar{U}) + \bar{\eth} + * (U)) J^2 + \eth (U) J \left(- \frac{1}{2 K} + \tfrac{1}{2} + * K\right)\right)\notag\\ + * &+ (1 - y) \bigg[\frac{1}{2} \left(- \frac{\partial_y (J)}{R} + \frac{2 + * \partial_{u} (R) \partial_y (\partial_y (J))}{R} + \partial_y + * (\partial_y (J)) W\right) + \partial_y (J) \left(\tfrac{1}{2} \partial_y (W) + * + \frac{1}{2 R}\right)\bigg]\notag\\ + * &+ (1 - y)^2 \bigg[ \frac{\partial_y (\partial_y (J)) }{4 R} \bigg], + * \f} + * where + * \f{align*} + * \mathcal{A}_J =& \tfrac{1}{4} \eth (\eth (\bar{J})) - \frac{1}{4 K^3} - + * \frac{\eth (\bar{\eth} (J \bar{J})) - (\eth (\bar{\eth} (\bar{J})) - 4 + \bar{J}) + * J}{16 K^3} + \frac{3}{4 K} - \frac{\eth (\bar{\eth} (\beta))}{4 K} \notag\\ + * &- \frac{\eth (\bar{\eth} (J)) \bar{J} (1 - \frac{1}{4 K^2})}{4 K} + + * \tfrac{1}{2} \eth (\bar{J}) \left(\eth (\beta) + \frac{\bar{\eth} (J \bar{J}) + * J}{4 K^3} - \frac{\bar{\eth} (J) (-1 + 2 K^2)}{4 K^3} - \tfrac{1}{2} + * Q\right)\\ + * \mathcal{B}_J =& - \frac{\eth (U) \bar{J} \partial_y (J \bar{J})}{4 K} + + * \tfrac{1}{2} \partial_y (W) + \frac{1}{4 R} + \tfrac{1}{4} \bar{\eth} (J) + * \partial_y (\bar{J}) U - \frac{\bar{\eth} (J \bar{J}) \partial_y (J \bar{J}) + * U}{8 K^2} \notag\\&- \tfrac{1}{4} J \partial_y (\eth (\bar{J})) \bar{U} + + * \tfrac{1}{4} (\eth (J \partial_y (\bar{J})) + \frac{J \eth (R) + * \partial_y(\bar{J})}{R}) \bar{U} \\ + * &+ (1 - y) \bigg[ \frac{\mathcal{D}_J \partial_{u} (R) + * \partial_y (J)}{R} - \tfrac{1}{4} \partial_y (J) \partial_y (\bar{J}) W + + * \frac{(\partial_y (J \bar{J}))^2 W}{16 K^2} \bigg] \\ + * &+ (1 - y)^2 \bigg[ - \frac{\partial_y (J) \partial_y (\bar{J})}{8 R} + + * \frac{(\partial_y (J \bar{J}))^2}{32 K^2 R} \bigg]\\ + * \mathcal{C}_J =& \tfrac{1}{2} \bar{\eth} (J) K (\eth (\beta) - \tfrac{1}{2} + * Q)\\ \mathcal{D}_J =& \tfrac{1}{4} \left(-2 \partial_y (\bar{J}) + + * \frac{\bar{J} \partial_y (J \bar{J})}{K^2}\right) + * \f} + */ +template <> +struct ComputeBondiIntegrand> { + public: + using pre_swsh_derivative_tags = + tmpl::list>, Tags::Dy, + Tags::Dy, Tags::Exp2Beta, Tags::BondiJ, + Tags::BondiQ, Tags::BondiU, Tags::BondiW>; + using swsh_derivative_tags = tmpl::list< + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Derivative< + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Eth>, + Spectral::Swsh::Tags::Derivative< + ::Tags::Multiplies, + Spectral::Swsh::Tags::EthEthbar>, + Spectral::Swsh::Tags::Derivative< + ::Tags::Multiplies, + Spectral::Swsh::Tags::Eth>, + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Derivative< + ::Tags::Multiplies>, + Spectral::Swsh::Tags::Eth>, + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Ethbar>, + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Derivative< + ::Tags::Multiplies>, + Spectral::Swsh::Tags::Ethbar>, + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Derivative, + Spectral::Swsh::Tags::Derivative>; + using integration_independent_tags = + tmpl::list; + using temporary_tags = + tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>, + std::integral_constant>, + ::Tags::SpinWeighted<::Tags::TempScalar<1, ComplexDataVector>, + std::integral_constant>, + ::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>, + std::integral_constant>>; + + using return_tags = + tmpl::append>, + temporary_tags>; + using argument_tags = + tmpl::append; + + template + static void apply( + const gsl::not_null>*> + regular_integrand_for_h, + const gsl::not_null>*> + script_aj, + const gsl::not_null>*> + script_bj, + const gsl::not_null>*> + script_cj, + const Args&... args) noexcept { + apply_impl(make_not_null(&get(*regular_integrand_for_h)), + make_not_null(&get(*script_aj)), make_not_null(&get(*script_bj)), + make_not_null(&get(*script_cj)), get(args)...); + } + + private: + static void apply_impl( + gsl::not_null*> + regular_integrand_for_h, + gsl::not_null*> script_aj, + gsl::not_null*> script_bj, + gsl::not_null*> script_cj, + const SpinWeighted& dy_dy_j, + const SpinWeighted& dy_j, + const SpinWeighted& dy_w, + const SpinWeighted& exp_2_beta, + const SpinWeighted& j, + const SpinWeighted& q, + const SpinWeighted& u, + const SpinWeighted& w, + const SpinWeighted& eth_beta, + const SpinWeighted& eth_eth_beta, + const SpinWeighted& eth_ethbar_beta, + const SpinWeighted& eth_ethbar_j, + const SpinWeighted& eth_ethbar_j_jbar, + const SpinWeighted& eth_j_jbar, + const SpinWeighted& eth_q, + const SpinWeighted& eth_u, + const SpinWeighted& eth_ubar_dy_j, + const SpinWeighted& ethbar_dy_j, + const SpinWeighted& ethbar_ethbar_j, + const SpinWeighted& ethbar_j, + const SpinWeighted& ethbar_jbar_dy_j, + const SpinWeighted& ethbar_jbar_q_minus_2_eth_beta, + const SpinWeighted& ethbar_q, + const SpinWeighted& ethbar_u, + const SpinWeighted& du_r_divided_by_r, + const SpinWeighted& eth_r_divided_by_r, + const SpinWeighted& k, + const SpinWeighted& one_minus_y, + const SpinWeighted& r) noexcept; +}; + +/*! + * \brief Computes the linear factor which multiplies \f$H\f$ in the + * equation which determines the radial (y) dependence of the Bondi quantity + * \f$H\f$. + * + * \details The quantity \f$H \equiv \partial_u J\f$ (evaluated at constant y) + * is defined via the Bondi form of the metric: + * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2 + * 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] + * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical + * angular dyad \f$q^A\f$: + * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B.\f] + * See \cite Bishop1997ik \cite Handmer2014qha for full details. + * + * We write the equations of motion in the compactified coordinate \f$ y \equiv + * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the + * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of + * the worldtube. The equation which determines \f$H\f$ on a surface of constant + * \f$u\f$ given \f$J\f$,\f$\beta\f$, \f$Q\f$, \f$U\f$, and \f$W\f$ on the same + * surface is written as + * \f[(1 - y) \partial_y H + H + (1 - y) J (\mathcal{D}_J + * H + \bar{\mathcal{D}}_J \bar{H}) = A_J + (1 - y) B_J.\f] + * The quantity \f$1 +(1 - y) J \mathcal{D}_J\f$ is the linear factor + * for the non-conjugated \f$H\f$, and is computed from the equation: + * \f[\mathcal{D}_J = \frac{1}{4}(-2 \partial_y \bar{J} + \frac{\bar{J} + * \partial_y (J \bar{J})}{K^2})\f] + */ +template <> +struct ComputeBondiIntegrand> { + public: + using pre_swsh_derivative_tags = + tmpl::list, Tags::BondiJ>; + using swsh_derivative_tags = tmpl::list<>; + using integration_independent_tags = tmpl::list; + using temporary_tags = + tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>, + std::integral_constant>>; + + using return_tags = tmpl::append>, + temporary_tags>; + using argument_tags = + tmpl::append; + + template + static void apply( + const gsl::not_null>*> + linear_factor_for_h, + const gsl::not_null>*> + script_djbar, + const Args&... args) noexcept { + apply_impl(make_not_null(&get(*linear_factor_for_h)), + make_not_null(&get(*script_djbar)), get(args)...); + } + + private: + static void apply_impl( + gsl::not_null*> linear_factor_for_h, + gsl::not_null*> script_djbar, + const SpinWeighted& dy_j, + const SpinWeighted& j, + const SpinWeighted& one_minus_y) noexcept; +}; + +/*! + * \brief Computes the linear factor which multiplies \f$\bar{H}\f$ in the + * equation which determines the radial (y) dependence of the Bondi quantity + * \f$H\f$. + * + * \details The quantity \f$H \equiv \partial_u J\f$ (evaluated at constant y) + * is defined via the Bondi form of the metric: + * \f[ds^2 = - \left(e^{2 \beta} (1 + r W) - r^2 h_{AB} U^A U^B\right) du^2 - 2 + * 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] + * Additional quantities \f$J\f$ and \f$K\f$ are defined using a spherical + * angular dyad \f$q^A\f$: + * \f[ J \equiv h_{A B} q^A q^B, K \equiv h_{A B} q^A \bar{q}^B.\f] + * See \cite Bishop1997ik \cite Handmer2014qha for full details. + * + * We write the equations of motion in the compactified coordinate \f$ y \equiv + * 1 - 2 R/ r\f$, where \f$r(u, \theta, \phi)\f$ is the Bondi radius of the + * \f$y=\f$ constant surface and \f$R(u,\theta,\phi)\f$ is the Bondi radius of + * the worldtube. The equation which determines \f$H\f$ on a surface of constant + * \f$u\f$ given \f$J\f$,\f$\beta\f$, \f$Q\f$, \f$U\f$, and \f$W\f$ on the same + * surface is written as + * \f[(1 - y) \partial_y H + H + (1 - y) J (\mathcal{D}_J H + + * \bar{\mathcal{D}}_J \bar{H}) = A_J + (1 - y) B_J.\f] + * The quantity \f$ (1 - y) J \bar{\mathcal{D}}_J\f$ is the linear factor + * for the non-conjugated \f$H\f$, and is computed from the equation: + * \f[\mathcal{D}_J = \frac{1}{4}(-2 \partial_y \bar{J} + \frac{\bar{J} + * \partial_y (J \bar{J})}{K^2})\f] + */ +template <> +struct ComputeBondiIntegrand> { + public: + using pre_swsh_derivative_tags = + tmpl::list, Tags::BondiJ>; + using swsh_derivative_tags = tmpl::list<>; + using integration_independent_tags = tmpl::list; + using temporary_tags = + tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<0, ComplexDataVector>, + std::integral_constant>>; + + using return_tags = + tmpl::append>, + temporary_tags>; + using argument_tags = + tmpl::append; + + template + static void apply( + const gsl::not_null>*> + linear_factor_for_conjugate_h, + const gsl::not_null>*> + script_djbar, + const Args&... args) noexcept { + apply_impl(make_not_null(&get(*linear_factor_for_conjugate_h)), + make_not_null(&get(*script_djbar)), get(args)...); + } + + private: + static void apply_impl( + gsl::not_null*> + linear_factor_for_conjugate_h, + gsl::not_null*> script_djbar, + const SpinWeighted& dy_j, + const SpinWeighted& j, + const SpinWeighted& one_minus_y) noexcept; +}; +} // namespace Cce diff --git a/src/Evolution/Systems/Cce/Tags.hpp b/src/Evolution/Systems/Cce/Tags.hpp new file mode 100644 index 000000000000..00d89ddc664f --- /dev/null +++ b/src/Evolution/Systems/Cce/Tags.hpp @@ -0,0 +1,233 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include "DataStructures/DataBox/DataBoxTag.hpp" +#include "DataStructures/SpinWeighted.hpp" +#include "DataStructures/Tensor/Tensor.hpp" +#include "DataStructures/Tensor/TypeAliases.hpp" + +namespace Cce { + +/// Tags for Cauchy Characteristic Extraction routines +namespace Tags { + +// Bondi parameter tags + +/// Bondi parameter \f$\beta\f$ +struct BondiBeta : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "BondiBeta"; } +}; + +/// \brief Bondi parameter \f$H = \partial_u J\f$. +/// \note The notation in the literature is not consistent regarding this +/// quantity, or whether it is denoted by an \f$H\f$ at all. The SpECTRE CCE +/// module consistently uses it to describe the (retarded) partial time +/// derivative of \f$J\f$ at fixed compactified radius \f$y\f$ (to be contrasted +/// with the physical Bondi radius, which is not directly used for numerical +/// grids). +struct BondiH : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "H"; } +}; + +/// Bondi parameter \f$J\f$ +struct BondiJ : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "J"; } +}; + +/// Bondi parameter \f$\bar{J}\f$ +struct BondiJbar : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "Jbar"; } +}; + +/// Bondi parameter \f$K = \sqrt{1 + J \bar{J}}\f$ +struct BondiK : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "K"; } +}; + +/// Bondi parameter \f$Q\f$ +struct BondiQ : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "Q"; } +}; + +/// Bondi parameter \f$\bar{Q}\f$ +struct BondiQbar : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "Qbar"; } +}; + +/// Bondi parameter \f$U\f$ +struct BondiU : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "U"; } +}; + +/// Bondi parameter \f$\bar{U}\f$ +struct BondiUbar : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "Ubar"; } +}; + +/// Bondi parameter \f$W\f$ +struct BondiW : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "W"; } +}; + +/// The derivative with respect to the numerical coordinate \f$y = 1 - 2R/r\f$, +/// where \f$R(u, \theta, \phi)\f$ is Bondi radius of the worldtube. +template +struct Dy : db::PrefixTag, db::SimpleTag { + using type = Scalar>; + using tag = Tag; + static const size_t dimension_to_differentiate = 2; + static const size_t dim = 3; + static std::string name() noexcept { return "Dy(" + Tag::name() + ")"; } +}; + +// prefix tags associated with the integrands which are used as input to solvers +// for the CCE equations + +/// A prefix tag representing a quantity that will appear on the right-hand side +/// of an explicitly regular differential equation +template +struct Integrand : db::PrefixTag, db::SimpleTag { + using type = Scalar>; + using tag = Tag; + static std::string name() noexcept { + return "Integrand(" + Tag::name() + ")"; + } +}; + +/// A prefix tag representing the boundary data for a quantity on the extraction +/// surface. +template +struct BoundaryValue : db::PrefixTag, db::SimpleTag { + using type = Scalar>; + using tag = Tag; + static std::string name() noexcept { + return "BoundaryValue(" + Tag::name() + ")"; + } +}; + +/// A prefix tag representing the coefficient of a pole part of the right-hand +/// side of a singular differential equation +template +struct PoleOfIntegrand : db::PrefixTag, db::SimpleTag { + using type = Scalar>; + using tag = Tag; + static std::string name() noexcept { + return "PoleOfIntegrand(" + Tag::name() + ")"; + } +}; + +/// A prefix tag representing the regular part of the right-hand side of a +/// regular differential equation +template +struct RegularIntegrand : db::PrefixTag, db::SimpleTag { + using type = Scalar>; + using tag = Tag; + static std::string name() noexcept { + return "RegularIntegrand(" + Tag::name() + ")"; + } +}; + +/// A prefix tag representing a linear factor that acts on `Tag`. To determine +/// the spin weight, It is assumed that the linear factor plays the role of +/// \f$L\f$ in an equation of the form, +/// \f$ (y - 1) \partial_y H + L H + L^\prime \bar{H} = A + (1 - y) B \f$ +template +struct LinearFactor : db::PrefixTag, db::SimpleTag { + using type = Scalar>; + using tag = Tag; + static std::string name() noexcept { + return "LinearFactor(" + Tag::name() + ")"; + } +}; + +/// A prefix tag representing a linear factor that acts on `Tag`. To determine +/// the spin weight, it is assumed that the linear factor plays the role of +/// \f$L^\prime\f$ in an equation of the form, +/// \f$ (y - 1) \partial_y H + L H + L^\prime \bar{H} = A + (1 - y) B \f$ +template +struct LinearFactorForConjugate : db::PrefixTag, db::SimpleTag { + using type = + Scalar>; + using tag = Tag; + static std::string name() noexcept { + return "LinearFactorForConjugate(" + Tag::name() + ")"; + } +}; + +// Below are additional tags for values which are frequently used in CCE +// calculations, and therefore worth caching + +/// Coordinate value \f$(1 - y)\f$, which will be cached and sent to the +/// implementing functions +struct OneMinusY : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "OneMinusY"; } +}; + +/// A tag for the first time derivative of the worldtube parameter +/// \f$\partial_u R\f$, where \f$R(u, \theta, \phi)\f$ is Bondi +/// radius of the worldtube. +struct DuR : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "DuR"; } +}; + +/// The value \f$\partial_u R / R\f$, where \f$R(u, \theta, \phi)\f$ is Bondi +/// radius of the worldtube. +struct DuRDividedByR : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "DuRDividedByR"; } +}; + +/// The value \f$\eth R / R\f$, where \f$R(u, \theta, \phi)\f$ is Bondi +/// radius of the worldtube. +struct EthRDividedByR : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "EthRDividedByR"; } +}; + +/// The value \f$\eth \eth R / R\f$, where \f$R(u, \theta, \phi)\f$ is Bondi +/// radius of the worldtube. +struct EthEthRDividedByR : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "EthEthRDividedByR"; } +}; + +/// The value \f$\eth \bar{\eth} R / R\f$, where \f$R(u, \theta, \phi)\f$ is +/// Bondi radius of the worldtube. +struct EthEthbarRDividedByR : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "EthEthbarRDividedByR"; } +}; + +/// The value \f$\exp(2\beta)\f$. +struct Exp2Beta : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "Exp2Beta"; } +}; + +/// The value \f$ \bar{J} (Q - 2 \eth \beta ) \f$. +struct JbarQMinus2EthBeta : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "JbarQMinus2EthBeta"; } +}; + +/// The Bondi radius \f$R(u, \theta, \phi)\f$ is of the worldtube. +struct BondiR : db::SimpleTag { + using type = Scalar>; + static std::string name() noexcept { return "R"; } +}; +} // namespace Tags +} // namespace Cce diff --git a/tests/Unit/DataStructures/Test_Tags.cpp b/tests/Unit/DataStructures/Test_Tags.cpp index 70215651f4a4..fca5aed2477f 100644 --- a/tests/Unit/DataStructures/Test_Tags.cpp +++ b/tests/Unit/DataStructures/Test_Tags.cpp @@ -7,16 +7,24 @@ #include #include "DataStructures/DataBox/DataBoxTag.hpp" +#include "DataStructures/SpinWeighted.hpp" #include "DataStructures/Tags.hpp" #include "DataStructures/Tensor/Tensor.hpp" #include "Utilities/TypeTraits.hpp" +class ComplexDataVector; class DataVector; class ModalVector; // IWYU pragma: no_forward_declare Tensor +// IWYU pragma: no_forward_declare SpinWeighted namespace { +struct ComplexScalarTag : db::SimpleTag { + using type = Scalar; + static std::string name() noexcept { return "ComplexScalar"; } +}; + struct ScalarTag : db::SimpleTag { using type = Scalar; static std::string name() noexcept { return "Scalar"; } @@ -59,9 +67,38 @@ void test_modal_tag() noexcept { CHECK(Tags::Modal>::name() == "Modal(I<1>)"); CHECK(Tags::Modal>::name() == "Modal(I<3>)"); } + +void test_spin_weighted_tag() noexcept { + static_assert( + cpp17::is_same_v< + typename Tags::SpinWeighted>::type, + Scalar>>, + "Failed testing Tags::SpinWeighted"); + CHECK(Tags::SpinWeighted>::name() == + "SpinWeighted(ComplexScalar, -2)"); +} + + +void test_multiplies_tag() noexcept { + using test_multiplies_tag = Tags::Multiplies< + Tags::SpinWeighted>, + Tags::SpinWeighted>>; + static_assert( + cpp17::is_same_v, + Scalar>>, + "Failed testing Tags::Multiplies for Tags::SpinWeighted operands"); + CHECK(test_multiplies_tag::name() == + "Multiplies(SpinWeighted(ComplexScalar, 1), " + "SpinWeighted(ComplexScalar, -2))"); +} } // namespace + SPECTRE_TEST_CASE("Unit.DataStructures.Tags", "[Unit][DataStructures]") { test_mean_tag(); test_modal_tag(); + test_spin_weighted_tag(); + test_multiplies_tag(); } diff --git a/tests/Unit/Evolution/Systems/CMakeLists.txt b/tests/Unit/Evolution/Systems/CMakeLists.txt index 1138f0a84ac4..720ef03012ed 100644 --- a/tests/Unit/Evolution/Systems/CMakeLists.txt +++ b/tests/Unit/Evolution/Systems/CMakeLists.txt @@ -2,6 +2,7 @@ # See LICENSE.txt for details. add_subdirectory(Burgers) +add_subdirectory(Cce) add_subdirectory(CurvedScalarWave) add_subdirectory(GeneralizedHarmonic) add_subdirectory(GrMhd) diff --git a/tests/Unit/Evolution/Systems/Cce/CMakeLists.txt b/tests/Unit/Evolution/Systems/Cce/CMakeLists.txt new file mode 100644 index 000000000000..b3e6daca3d78 --- /dev/null +++ b/tests/Unit/Evolution/Systems/Cce/CMakeLists.txt @@ -0,0 +1,15 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +set(LIBRARY "Test_Cce") + +set(LIBRARY_SOURCES + Test_Equations.cpp + ) + +add_test_library( + ${LIBRARY} + "Evolution/Systems/Cce/" + "${LIBRARY_SOURCES}" + "Cce" + ) diff --git a/tests/Unit/Evolution/Systems/Cce/Equations.py b/tests/Unit/Evolution/Systems/Cce/Equations.py new file mode 100644 index 000000000000..c26c458c2f12 --- /dev/null +++ b/tests/Unit/Evolution/Systems/Cce/Equations.py @@ -0,0 +1,205 @@ +# Distributed under the MIT License. +# See LICENSE.txt for details. + +import numpy as np + +# Test function for beta equation + + +def integrand_for_beta(dy_j, j, one_minus_y): + dy_jbar = np.conj(dy_j) + jbar = np.conj(j) + dy_j_jbar = j * dy_jbar + jbar * dy_j + k_squared = 1.0 + j * jbar + return one_minus_y / 8.0 * (dy_j * dy_jbar \ + - dy_j_jbar**2 / (4.0 * k_squared)) + + +# Test functions for Q equation + + +def integrand_for_q_pole_part(eth_beta): + return -4.0 * eth_beta + + +def integrand_for_q_regular_part( + _, dy_beta, dy_j, j, eth_dy_beta, eth_j_jbar, eth_jbar_dy_j, + ethbar_dy_j, ethbar_j, eth_r_divided_by_r, k): + # script_aq input is unused, needed only to match function signature from + # the C++ + eth_dy_jbar = np.conj(ethbar_dy_j) + dy_jbar = np.conj(dy_j) + ethbar_r_divided_by_r = np.conj(eth_r_divided_by_r) + eth_jbar = np.conj(ethbar_j) + jbar = np.conj(j) + script_aq = - eth_jbar_dy_j / 4.0 + j * eth_dy_jbar / 4.0 \ + - eth_jbar * dy_j / 4.0 \ + + eth_j_jbar * (jbar * dy_j + j * dy_jbar)\ + / (8.0 * (1.0 + j * jbar)) \ + + (j * dy_jbar - jbar * dy_j) * eth_r_divided_by_r / 4.0 + return - (2.0 * script_aq + 2.0 * j * np.conj(script_aq) / k \ + - 2.0 * eth_dy_beta + ethbar_dy_j / k \ + - 2.0 * dy_beta * eth_r_divided_by_r \ + + dy_j * ethbar_r_divided_by_r / k) + + +# Test function for U equation + + +def integrand_for_u(exp_2_beta, j, q, k, r): + qbar = np.conj(q) + return exp_2_beta / (2.0 * r) * (k * q - j * qbar) + + +# Test functions for W equation + + +def integrand_for_w_pole_part(ethbar_u): + eth_ubar = np.conj(ethbar_u) + return eth_ubar + ethbar_u + + +def integrand_for_w_regular_part( + _, dy_u, exp_2_beta, j, q, eth_beta, eth_eth_beta, + eth_ethbar_beta, eth_ethbar_j, eth_ethbar_j_jbar, eth_j_jbar, + ethbar_dy_u, ethbar_ethbar_j, ethbar_j, eth_r_divided_by_r, k, r): + # script_av input is unused, needed only to match function signature from + # the C++ + dy_ubar = np.conj(dy_u) + ethbar_beta = np.conj(eth_beta) + eth_dy_ubar = np.conj(ethbar_dy_u) + ethbar_ethbar_beta = np.conj(eth_eth_beta) + eth_eth_jbar = np.conj(ethbar_ethbar_j) + eth_jbar = np.conj(ethbar_j) + ethbar_j_jbar = np.conj(eth_j_jbar) + ethbar_r_divided_by_r = np.conj(eth_r_divided_by_r) + jbar = np.conj(j) + k_squared = (1.0 + j * jbar) + qbar = np.conj(q) + script_av = eth_beta * eth_jbar + ethbar_ethbar_j / 2.0 \ + + j * ethbar_beta**2 + j * ethbar_ethbar_beta \ + + eth_j_jbar * ethbar_j_jbar / (8.0 * k_squared * k) \ + + 1 / (2.0 * k) - eth_ethbar_j_jbar / (8.0 * k) \ + - eth_j_jbar * ethbar_beta / (2.0 * k) \ + - eth_jbar * ethbar_j / (4.0 * k) \ + - eth_ethbar_j * jbar / (4.0 * k) + k / 2.0 \ + - eth_ethbar_beta * k - eth_beta * ethbar_beta * k \ + + 1.0 / 4.0 * (- k * q * qbar + j * qbar**2) + return 1.0 / 4.0 * eth_dy_ubar + 1.0 / 4.0 * ethbar_dy_u \ + + 1.0 / 4.0 * dy_u * ethbar_r_divided_by_r \ + + 1.0 / 4.0 * dy_ubar * eth_r_divided_by_r - 1.0 / (2.0 * r) \ + + exp_2_beta * (script_av + np.conj(script_av)) / (4.0 * r) + + +# Test functions for H equation + + +def integrand_for_h_pole_part(j, u, w, eth_u, ethbar_j, ethbar_jbar_u, + ethbar_u, k): + eth_ubar = np.conj(ethbar_u) + eth_j_ubar = np.conj(ethbar_jbar_u) + return - 1.0 / 2.0 * eth_j_ubar - j * eth_ubar - 1.0 / 2.0 * j * ethbar_u \ + - k * eth_u - 1.0 / 2.0 * u * ethbar_j + 2.0 * j * w + + +def integrand_for_h_regular_part( + _0, _1, _2, dy_dy_j, dy_j, dy_w, exp_2_beta, j, q, + u, w, eth_beta, eth_eth_beta, eth_ethbar_beta, eth_ethbar_j, + eth_ethbar_j_jbar, eth_j_jbar, eth_q, eth_u, eth_ubar_dy_j, + ethbar_dy_j, ethbar_ethbar_j, ethbar_j, ethbar_jbar_dy_j, + ethbar_jbar_q_minus_2_eth_beta, ethbar_q, ethbar_u, du_r_divided_by_r, + eth_r_divided_by_r, k, one_minus_y, r): + # script_aj, script_bj, and script_cj input is unused, needed only to match + # function signature from the C++ + jbar = np.conj(j) + k_squared = 1.0 + j * jbar + eth_eth_jbar = np.conj(ethbar_ethbar_j) + ethbar_eth_jbar = np.conj(eth_ethbar_j) + eth_jbar = np.conj(ethbar_j) + ethbar_j_jbar = np.conj(eth_j_jbar) + + script_aj = 1.0 / 4.0 * eth_eth_jbar - 1.0 / (4.0 * k * k_squared) \ + - eth_ethbar_j_jbar / (16.0 * k * k_squared) \ + + j * ethbar_eth_jbar / (16.0 * k * k_squared) \ + + 3.0 / (4.0 * k) - eth_ethbar_beta / (4.0 * k) \ + - eth_ethbar_j * jbar * (1.0 - 1.0 / (4.0 * k_squared)) \ + / (4.0 * k) \ + + 1.0 / 2.0 * eth_jbar * (eth_beta + j * ethbar_j_jbar \ + / (4.0 * k * k_squared) \ + - ethbar_j * (-1.0 + 2.0 * k_squared) \ + / (4.0 * k * k_squared) - 1.0 / 2.0 * q) + + dy_jbar = np.conj(dy_j) + dy_j_jbar = j * dy_jbar + jbar * dy_j + eth_dy_jbar = np.conj(ethbar_dy_j) + eth_j_dy_jbar = np.conj(ethbar_jbar_dy_j) + ubar = np.conj(u) + + script_bj = (- eth_u * jbar * dy_j_jbar / (4.0 * k) + 1.0 / 2.0 * dy_w \ + + 1.0 / (4.0 * r) + 1.0 / 4.0 * ethbar_j * dy_jbar * u \ + - ethbar_j_jbar * dy_j_jbar * u / (8.0 * k_squared) \ + - 1.0 / 4.0 * j * eth_dy_jbar * ubar \ + + 1.0 / 4.0 * eth_j_dy_jbar * ubar) \ + + one_minus_y \ + * (du_r_divided_by_r * dy_j * 1.0 / 4.0 \ + * (-2.0 * dy_jbar + jbar * dy_j_jbar / k_squared) \ + - 1.0/4.0 * dy_j * dy_jbar * w \ + + w * dy_j_jbar**2 / (16.0 * k_squared)) \ + + one_minus_y**2 * (- dy_j * dy_jbar / (8.0 * r) \ + + dy_j_jbar**2 / (32.0 * k_squared * r)) + + script_cj = 1.0 / 2.0 * ethbar_j * k * (eth_beta - 1.0 / 2.0 * q) + + eth_j_qbar_minus_2_ethbar_beta = np.conj(ethbar_jbar_q_minus_2_eth_beta) + eth_qbar = np.conj(ethbar_q) + ethbar_r_divided_by_r = np.conj(eth_r_divided_by_r) + ethbar_ubar = np.conj(eth_u) + eth_ubar = np.conj(ethbar_u) + ubar = np.conj(u) + dy_jbar = np.conj(dy_j) + + return -1.0 / 2.0 * eth_ubar_dy_j - 1.0 / 2.0 * u * ethbar_dy_j \ + - 1.0 / 2.0 * u * dy_j * ethbar_r_divided_by_r \ + + j * (script_bj + np.conj(script_bj)) \ + + exp_2_beta / (2.0 * r) \ + * (script_cj + j**2 / k_squared * np.conj(script_cj) \ + + eth_eth_beta - 1.0 / 2.0 * eth_q \ + - j * (script_aj + np.conj(script_aj)) \ + + eth_j_qbar_minus_2_ethbar_beta / (4.0 * k) \ + - j * eth_qbar / (4.0 * k) + (eth_beta - 1.0 / 2.0 * q)**2) \ + - dy_j \ + * (jbar * eth_u / (2.0 * k) - 1.0 / 2.0 * ethbar_ubar * j * k \ + + 1.0 / 4.0 * (eth_ubar - ethbar_u) * k_squared + \ + + 1.0 / 2.0 * eth_r_divided_by_r * ubar - 1.0 / 2.0 * w) \ + + dy_jbar \ + * (-1.0 / 4.0 * j**2 * (- eth_ubar + ethbar_u) \ + + eth_u * j * (j * jbar /( 2.0 * k))) \ + + one_minus_y \ + * (1.0 / 2.0 * (-dy_j / r + 2.0 * du_r_divided_by_r * dy_dy_j \ + + w * dy_dy_j) + dy_j * (1.0 / 2.0 * dy_w \ + + 1.0 / (2.0 * r))) \ + + one_minus_y**2 * dy_dy_j / (4.0 * r) + + +def linear_factor_for_h(_, dy_j, j, one_minus_y): + # script_dj input is unused, needed only to match function signature from + # the C++ + dy_jbar = np.conj(dy_j) + jbar = np.conj(j) + k_squared = 1.0 + j * np.conj(j) + + return 1.0 + ((1.0 / 4.0) * np.conj(one_minus_y) * j \ + * (-2.0 * dy_jbar + jbar * (jbar * dy_j + + j * dy_jbar) / k_squared)) + + +def linear_factor_for_conjugate_h(_, dy_j, j, one_minus_y): + # script_dj input is unused, needed only to match function signature from + # the C++ + dy_jbar = np.conj(dy_j) + jbar = np.conj(j) + k_squared = 1.0 + j * np.conj(j) + + return (1.0 / 4.0) * one_minus_y * j \ + * np.conj(-2.0 * dy_jbar + jbar \ + * (jbar * dy_j + j * dy_jbar) / k_squared) diff --git a/tests/Unit/Evolution/Systems/Cce/Test_Equations.cpp b/tests/Unit/Evolution/Systems/Cce/Test_Equations.cpp new file mode 100644 index 000000000000..e4ec056a269e --- /dev/null +++ b/tests/Unit/Evolution/Systems/Cce/Test_Equations.cpp @@ -0,0 +1,182 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "tests/Unit/TestingFramework.hpp" + +#include +#include +#include + +#include "DataStructures/ComplexDataVector.hpp" +#include "DataStructures/DataBox/DataBox.hpp" +#include "DataStructures/DataBox/DataBoxTag.hpp" +#include "DataStructures/Variables.hpp" // IWYU pragma: keep +#include "Evolution/Systems/Cce/Equations.hpp" +#include "Evolution/Systems/Cce/Tags.hpp" // IWYU pragma: keep +#include "Utilities/Gsl.hpp" +#include "Utilities/TMPL.hpp" +#include "tests/Unit/Pypp/CheckWithRandomValues.hpp" +#include "tests/Unit/Pypp/SetupLocalPythonEnvironment.hpp" + +// IWYU pragma: no_forward_declare Cce::Tags::BondiBeta +// IWYU pragma: no_forward_declare Cce::Tags::H +// IWYU pragma: no_forward_declare Cce::Tags::Q +// IWYU pragma: no_forward_declare Cce::Tags::U +// IWYU pragma: no_forward_declare Cce::Tags::W +// IWYU pragma: no_forward_declare Cce::Tags::Integrand +// IWYU pragma: no_forward_declare Cce::Tags::LinearFactor +// IWYU pragma: no_forward_declare Cce::Tags::LinearFactorForConjugate +// IWYU pragma: no_forward_declare Cce::Tags::PoleOfIntegrand +// IWYU pragma: no_forward_declare Cce::Tags::RegularIntegrand +// IWYU pragma: no_forward_declare Cce::TagsCategory::IndependentOfBondiIntegration +// IWYU pragma: no_forward_declare Cce::TagsCategory::SwshDerivatives +// IWYU pragma: no_forward_declare Cce::TagsCategory::TagsToSwshDifferentiate +// IWYU pragma: no_forward_declare Cce::TagsCategory::Temporary +// IWYU pragma: no_forward_declare Cce::ComputeBondiIntegrand +// IWYU pragma: no_forward_declare Variables + +namespace Cce { + +namespace { + +template +using ForwardIgnoreIndex = Type; + +// This is a wrapper struct such that the evaluate function can be called with +// the full set of arguments that would ordinarily be packed into the +// DataBox. This is necessary for constructing a function that can be compared +// to the pypp versions. +template +struct MutationFromArguments; + +template +struct MutationFromArguments, + NumberOfGridPoints, Args...> { + static ComplexDataVector evaluate(const Args&... args) noexcept { + auto box = db::create>( + db::item_type{typename db::item_type::type{ + ComplexDataVector{NumberOfGridPoints}}}, + db::item_type{ + typename db::item_type::type{args}}...); + db::mutate_apply(make_not_null(&box)); + return ComplexDataVector{get(db::get(box)).data()}; + } +}; + +// Creates a proxy struct which acts like the function being tested, but with +// arguments not packed into a DataBox, then compares the result of that +// function to the result of a python function +template +void forward_to_pypp_with(std::index_sequence /*meta*/, + const std::string& python_function) noexcept { + using Evaluatable = + MutationFromArguments...>; + + pypp::check_with_random_values<1>(&(Evaluatable::evaluate), "Equations", + python_function, {{{-10.0, 10.0}}}, + DataType{NumberOfGridPoints}); +} + +template +using all_arguments_for_integrand = + tmpl::append::temporary_tags, + typename ComputeBondiIntegrand::argument_tags>; + +template +struct python_function_for_bondi_integrand; + +template <> +struct python_function_for_bondi_integrand> { + static std::string name() noexcept { return "integrand_for_beta"; } +}; +template <> +struct python_function_for_bondi_integrand< + Tags::PoleOfIntegrand> { + static std::string name() noexcept { + return "integrand_for_q_pole_part"; + } +}; +template <> +struct python_function_for_bondi_integrand< + Tags::RegularIntegrand> { + static std::string name() noexcept { + return "integrand_for_q_regular_part"; + } +}; +template <> +struct python_function_for_bondi_integrand> { + static std::string name() noexcept { return "integrand_for_u"; } +}; +template <> +struct python_function_for_bondi_integrand< + Tags::PoleOfIntegrand> { + static std::string name() noexcept { + return "integrand_for_w_pole_part"; + } +}; +template <> +struct python_function_for_bondi_integrand< + Tags::RegularIntegrand> { + static std::string name() noexcept { + return "integrand_for_w_regular_part"; + } +}; +template <> +struct python_function_for_bondi_integrand< + Tags::PoleOfIntegrand> { + static std::string name() noexcept { + return "integrand_for_h_pole_part"; + } +}; +template <> +struct python_function_for_bondi_integrand< + Tags::RegularIntegrand> { + static std::string name() noexcept { + return "integrand_for_h_regular_part"; + } +}; +template <> +struct python_function_for_bondi_integrand> { + static std::string name() noexcept { return "linear_factor_for_h"; } +}; +template <> +struct python_function_for_bondi_integrand< + Tags::LinearFactorForConjugate> { + static std::string name() noexcept { + return "linear_factor_for_conjugate_h"; + } +}; + +SPECTRE_TEST_CASE("Unit.Evolution.Systems.Cce.Equations", "[Unit][Evolution]") { + pypp::SetupLocalPythonEnvironment local_python_env{"Evolution/Systems/Cce/"}; + + using all_bondi_tags = tmpl::list; + tmpl::for_each([](auto x) { + using bondi_tag = typename decltype(x)::type; + tmpl::for_each>( + [](auto y) { + using bondi_integrand_tag = typename decltype(y)::type; + using bondi_integrand_argument_list = + all_arguments_for_integrand; + // We check the equations with the DataBox interface, as it is easy + // to generate the box with the appropriate type lists, + // which are then used to create wrappers that can be compared to Pypp + // equations without explicitly enumerating the long argument lists in + // this test. + forward_to_pypp_with, + bondi_integrand_argument_list, + bondi_integrand_tag, 5, ComplexDataVector>( + std::make_index_sequence< + tmpl::size::value>{}, + python_function_for_bondi_integrand::name()); + }); + }); +} +} // namespace +} // namespace Cce