Enable definite_integral for spherical harmonic basis

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<double>& 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];
+        }
       }
     }
   }

tests/Unit/Helpers/NumericalAlgorithms/SphericalHarmonics/CMakeLists.txt

@@ -14,6 +14,7 @@ target_link_libraries(
   ${LIBRARY}
   PUBLIC
   DataStructures
+  Spectral
   SphericalHarmonics
   Utilities
   )

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<double>(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<double>(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<double>(pow_nx_ + pow_ny_) *
+              pow(sin(theta_), pow_nx_ + pow_ny_ - 1) *
+              pow(cos(theta_), pow_nz_ + 1) -
+          static_cast<double>(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<double>(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<double>(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<double>(pow_ny_) * pow(cos(phi_), pow_nx_ + 1) *
+              pow(sin(phi_), pow_ny_ - 1) -
+          static_cast<double>(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<double>(falling_factorial(pow_nx_, pow_nx_ / 2)) *
+         static_cast<double>(falling_factorial(pow_ny_, pow_ny_ / 2)) *
+         static_cast<double>(falling_factorial(pow_nz_, pow_nz_ / 2)) /
+         static_cast<double>(
+             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)};

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 <size_t l, int m>
 class Ylm {
  public:

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<Spectral::Basis::FiniteDifference,