diff --git a/src/NumericalAlgorithms/LinearOperators/DefiniteIntegral.cpp b/src/NumericalAlgorithms/LinearOperators/DefiniteIntegral.cpp index 2624bbe228d0..8a5a18cc304a 100644 --- a/src/NumericalAlgorithms/LinearOperators/DefiniteIntegral.cpp +++ b/src/NumericalAlgorithms/LinearOperators/DefiniteIntegral.cpp @@ -9,6 +9,8 @@ #include "DataStructures/DataVector.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" #include "NumericalAlgorithms/Spectral/Spectral.hpp" +#include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp" +#include "NumericalAlgorithms/SphericalHarmonics/SpherepackCache.hpp" #include "Utilities/Blas.hpp" #include "Utilities/ErrorHandling/Assert.hpp" @@ -82,33 +84,51 @@ double definite_integral<3>(const DataVector& integrand, const Mesh<3>& mesh) { const auto sliced_meshes = mesh.slices(); const size_t x_size = sliced_meshes[0].number_of_grid_points(); const size_t x_last_unrolled = x_size - x_size % 2; - const size_t y_size = sliced_meshes[1].number_of_grid_points(); - const size_t z_size = sliced_meshes[2].number_of_grid_points(); const double* const w_x = Spectral::quadrature_weights(sliced_meshes[0]).data(); - const double* const w_y = - Spectral::quadrature_weights(sliced_meshes[1]).data(); - const double* const w_z = - Spectral::quadrature_weights(sliced_meshes[2]).data(); - double result = 0.0; - for (size_t k = 0; k < z_size; ++k) { - for (size_t j = 0; j < y_size; ++j) { - // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) - const double prod = w_z[k] * w_y[j]; - const size_t offset = x_size * (j + y_size * k); - // Unrolling at 2 gives better performance on AVX machines (the most - // common in 2018). The stride will probably need to be updated as - // hardware changes. Note: using a single loop is ~15% faster when - // x_size == 3, and both loop styles are comparable for x_size == 2. + if (mesh.basis(1) == Spectral::Basis::SphericalHarmonic) { + const auto& ylm = ylm::get_spherepack_cache(mesh.extents(1) - 1); + const size_t angular_size = ylm.physical_size(); + const std::vector& w_ylm = ylm.integration_weights(); + for (size_t j = 0; j < angular_size; ++j) { + const size_t offset = j * x_size; for (size_t i = 0; i < x_last_unrolled; i += 2) { // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) - result += prod * w_x[i] * integrand[i + offset] + - prod * w_x[i + 1] * integrand[i + 1 + offset]; // NOLINT + result += w_ylm[j] * w_x[i] * integrand[i + offset] + + w_ylm[j] * w_x[i + 1] * integrand[i + 1 + offset]; // NOLINT } for (size_t i = x_last_unrolled; i < x_size; ++i) { // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) - result += prod * w_x[i] * integrand[i + offset]; + result += w_ylm[j] * w_x[i] * integrand[i + offset]; + } + } + } else { + const size_t y_size = sliced_meshes[1].number_of_grid_points(); + const size_t z_size = sliced_meshes[2].number_of_grid_points(); + const double* const w_y = + Spectral::quadrature_weights(sliced_meshes[1]).data(); + const double* const w_z = + Spectral::quadrature_weights(sliced_meshes[2]).data(); + + for (size_t k = 0; k < z_size; ++k) { + for (size_t j = 0; j < y_size; ++j) { + // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) + const double prod = w_z[k] * w_y[j]; + const size_t offset = x_size * (j + y_size * k); + // Unrolling at 2 gives better performance on AVX machines (the most + // common in 2018). The stride will probably need to be updated as + // hardware changes. Note: using a single loop is ~15% faster when + // x_size == 3, and both loop styles are comparable for x_size == 2. + for (size_t i = 0; i < x_last_unrolled; i += 2) { + // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) + result += prod * w_x[i] * integrand[i + offset] + + prod * w_x[i + 1] * integrand[i + 1 + offset]; // NOLINT + } + for (size_t i = x_last_unrolled; i < x_size; ++i) { + // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic) + result += prod * w_x[i] * integrand[i + offset]; + } } } } diff --git a/tests/Unit/Helpers/NumericalAlgorithms/SphericalHarmonics/CMakeLists.txt b/tests/Unit/Helpers/NumericalAlgorithms/SphericalHarmonics/CMakeLists.txt index a66a0f663f42..9883aac190fc 100644 --- a/tests/Unit/Helpers/NumericalAlgorithms/SphericalHarmonics/CMakeLists.txt +++ b/tests/Unit/Helpers/NumericalAlgorithms/SphericalHarmonics/CMakeLists.txt @@ -14,6 +14,7 @@ target_link_libraries( ${LIBRARY} PUBLIC DataStructures + Spectral SphericalHarmonics Utilities ) diff --git a/tests/Unit/Helpers/NumericalAlgorithms/SphericalHarmonics/YlmTestFunctions.cpp b/tests/Unit/Helpers/NumericalAlgorithms/SphericalHarmonics/YlmTestFunctions.cpp index d8ac2478b497..80f0be1bb499 100644 --- a/tests/Unit/Helpers/NumericalAlgorithms/SphericalHarmonics/YlmTestFunctions.cpp +++ b/tests/Unit/Helpers/NumericalAlgorithms/SphericalHarmonics/YlmTestFunctions.cpp @@ -8,6 +8,87 @@ namespace YlmTestFunctions { +ProductOfPolynomials::ProductOfPolynomials(const size_t n_r, const size_t L, + const size_t M, const size_t pow_nx, + const size_t pow_ny, + const size_t pow_nz) + : pow_nx_{pow_nx}, pow_ny_{pow_ny}, pow_nz_{pow_nz} { + ASSERT(std::hypot(pow_nx_, pow_ny_, pow_nz_) <= L, + "Cannot represent function on mesh"); + ASSERT(std::hypot(pow_nx_, pow_ny_) <= M, + "Cannot represent function on mesh"); + const Mesh<3> mesh{ + {n_r, L + 1, 2 * M + 1}, + {Spectral::Basis::Legendre, Spectral::Basis::SphericalHarmonic, + Spectral::Basis::SphericalHarmonic}, + {Spectral::Quadrature::Gauss, Spectral::Quadrature::Gauss, + Spectral::Quadrature::Equiangular}}; + const auto xi = logical_coordinates(mesh); + n_pts_ = mesh.number_of_grid_points(); + theta_ = xi.get(1); + phi_ = xi.get(2); +} + +DataVector ProductOfPolynomials::f() const { + return pow(sin(theta_), pow_nx_ + pow_ny_) * pow(cos(theta_), pow_nz_) * + pow(cos(phi_), pow_nx_) * pow(sin(phi_), pow_ny_); +} + +DataVector ProductOfPolynomials::df_dth() const { + if (pow_nx_ + pow_ny_ + pow_nz_ == 0) { + return DataVector{n_pts_, 0.0}; + } + if (pow_nx_ + pow_ny_ == 0) { + return -static_cast(pow_nz_) * sin(theta_) * + pow(cos(theta_), pow_nz_ - 1) * pow(cos(phi_), pow_nx_) * + pow(sin(phi_), pow_ny_); + } + if (pow_nz_ == 0) { + return static_cast(pow_nx_ + pow_ny_) * + pow(sin(theta_), pow_nx_ + pow_ny_ - 1) * cos(theta_) * + pow(cos(phi_), pow_nx_) * pow(sin(phi_), pow_ny_); + } + return (static_cast(pow_nx_ + pow_ny_) * + pow(sin(theta_), pow_nx_ + pow_ny_ - 1) * + pow(cos(theta_), pow_nz_ + 1) - + static_cast(pow_nz_) * + pow(sin(theta_), pow_nx_ + pow_ny_ + 1) * + pow(cos(theta_), pow_nz_ - 1)) * + pow(cos(phi_), pow_nx_) * pow(sin(phi_), pow_ny_); +} + +DataVector ProductOfPolynomials::df_dph() const { + if (pow_nx_ + pow_ny_ == 0) { + return DataVector{n_pts_, 0.0}; + } + if (pow_nx_ == 0) { + return static_cast(pow_ny_) * pow(sin(theta_), pow_nx_ + pow_ny_) * + pow(cos(theta_), pow_nz_) * cos(phi_) * pow(sin(phi_), pow_ny_ - 1); + } + if (pow_ny_ == 0) { + return -static_cast(pow_nx_) * pow(sin(theta_), pow_nx_ + pow_ny_) * + pow(cos(theta_), pow_nz_) * pow(cos(phi_), pow_nx_ - 1) * sin(phi_); + } + return pow(sin(theta_), pow_nx_ + pow_ny_) * pow(cos(theta_), pow_nz_) * + (static_cast(pow_ny_) * pow(cos(phi_), pow_nx_ + 1) * + pow(sin(phi_), pow_ny_ - 1) - + static_cast(pow_nx_) * *pow(cos(phi_), pow_nx_ - 1) * + pow(sin(phi_), pow_ny_ + 1)); +} + +double ProductOfPolynomials::definite_integral() const { + if ((pow_nx_ % 2 == 1) or (pow_ny_ % 2 == 1) or (pow_nz_ % 2 == 1)) { + return 0.0; + } + return 4.0 * M_PI * + static_cast(falling_factorial(pow_nx_, pow_nx_ / 2)) * + static_cast(falling_factorial(pow_ny_, pow_ny_ / 2)) * + static_cast(falling_factorial(pow_nz_, pow_nz_ / 2)) / + static_cast( + falling_factorial(pow_nx_ + pow_ny_ + pow_nz_ + 1, + (pow_nx_ + pow_ny_ + pow_nz_) / 2 + 1)); +} + template <> DataVector Ylm<0, 0>::f() const { return DataVector{n_pts_, 1.0 / sqrt(4.0 * M_PI)}; diff --git a/tests/Unit/Helpers/NumericalAlgorithms/SphericalHarmonics/YlmTestFunctions.hpp b/tests/Unit/Helpers/NumericalAlgorithms/SphericalHarmonics/YlmTestFunctions.hpp index 1aded0e46596..0b313bfe758f 100644 --- a/tests/Unit/Helpers/NumericalAlgorithms/SphericalHarmonics/YlmTestFunctions.hpp +++ b/tests/Unit/Helpers/NumericalAlgorithms/SphericalHarmonics/YlmTestFunctions.hpp @@ -18,6 +18,25 @@ namespace YlmTestFunctions { +class ProductOfPolynomials { + public: + ProductOfPolynomials(size_t n_r, size_t L, size_t M, size_t pow_nx, + size_t pow_ny, size_t pow_nz); + DataVector f() const; + DataVector df_dth() const; + // This is the Pfaffiaan derivative (extra factor 1/sin(th)) + DataVector df_dph() const; + double definite_integral() const; + + private: + size_t n_pts_; + size_t pow_nx_; + size_t pow_ny_; + size_t pow_nz_; + DataVector theta_; + DataVector phi_; +}; + template class Ylm { public: diff --git a/tests/Unit/NumericalAlgorithms/LinearOperators/Test_DefiniteIntegral.cpp b/tests/Unit/NumericalAlgorithms/LinearOperators/Test_DefiniteIntegral.cpp index 484588c544bb..cffcdba318ee 100644 --- a/tests/Unit/NumericalAlgorithms/LinearOperators/Test_DefiniteIntegral.cpp +++ b/tests/Unit/NumericalAlgorithms/LinearOperators/Test_DefiniteIntegral.cpp @@ -9,8 +9,10 @@ #include "DataStructures/Index.hpp" #include "DataStructures/IndexIterator.hpp" #include "Framework/TestHelpers.hpp" +#include "Helpers/NumericalAlgorithms/SphericalHarmonics/YlmTestFunctions.hpp" #include "NumericalAlgorithms/LinearOperators/DefiniteIntegral.hpp" #include "NumericalAlgorithms/Spectral/Basis.hpp" +#include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp" #include "NumericalAlgorithms/Spectral/Mesh.hpp" #include "NumericalAlgorithms/Spectral/Quadrature.hpp" #include "NumericalAlgorithms/Spectral/Spectral.hpp" @@ -105,6 +107,35 @@ void test_definite_integral_3d(const Mesh<3>& mesh) { } } } + +void test_definite_integral_spherical_shell(const size_t n_r, const size_t L) { + CAPTURE(n_r); + CAPTURE(L); + const Mesh<3> mesh{ + {n_r, L + 1, 2 * L + 1}, + {Spectral::Basis::Legendre, Spectral::Basis::SphericalHarmonic, + Spectral::Basis::SphericalHarmonic}, + {Spectral::Quadrature::Gauss, Spectral::Quadrature::Gauss, + Spectral::Quadrature::Equiangular}}; + const auto xi_vector = logical_coordinates(mesh); + const DataVector r = xi_vector.get(0) + 2.0; + for (size_t pow_nx = 0; pow_nx <= L; ++pow_nx) { + CAPTURE(pow_nx); + for (size_t pow_ny = 0; pow_ny <= L - pow_nx; ++pow_ny) { + CAPTURE(pow_ny); + for (size_t pow_nz = 0; pow_nz <= L - pow_nx - pow_ny; ++pow_nz) { + CAPTURE(pow_nz); + const YlmTestFunctions::ProductOfPolynomials y_lm{ + n_r, L, L, pow_nx, pow_ny, pow_nz}; + + const DataVector integrand = r * y_lm.f(); + CHECK(4.0 * y_lm.definite_integral() == + approx(definite_integral(integrand, mesh))); + } + } + } +} + } // namespace SPECTRE_TEST_CASE("Unit.Numerical.LinearOperators.DefiniteIntegral", @@ -137,6 +168,12 @@ SPECTRE_TEST_CASE("Unit.Numerical.LinearOperators.DefiniteIntegral", } } + for (size_t n_r = 2; n_r < 5; ++n_r) { + for (size_t L = 2; L < 9; ++L) { + test_definite_integral_spherical_shell(n_r, L); + } + } + // Test finite difference integral constexpr size_t min_extents_fd = Spectral::minimum_number_of_points