Support logical_coordinates for SphericalHarmonic Basis

src/NumericalAlgorithms/Spectral/CMakeLists.txt

@@ -40,10 +40,12 @@ target_link_libraries(
   DataStructures
   ErrorHandling
   Options
+  Utilities
   PRIVATE
   BLAS::BLAS
   LAPACK::LAPACK
   RootFinding
+  SPHEREPACK
   )
 
 add_subdirectory(Python)

src/NumericalAlgorithms/Spectral/LogicalCoordinates.cpp

@@ -14,6 +14,7 @@
 #include "NumericalAlgorithms/Spectral/Spectral.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/SetNumberOfGridPoints.hpp"
+#include "Utilities/Spherepack.hpp"
 
 template <size_t VolumeDim>
 void logical_coordinates(
@@ -22,11 +23,64 @@ void logical_coordinates(
     const Mesh<VolumeDim>& mesh) {
   set_number_of_grid_points(logical_coords, mesh.number_of_grid_points());
   for (size_t d = 0; d < VolumeDim; ++d) {
-    const auto& collocation_points_in_this_dim =
-        Spectral::collocation_points(mesh.slice_through(d));
-    for (IndexIterator<VolumeDim> index(mesh.extents()); index; ++index) {
-      logical_coords->get(d)[index.collapsed_index()] =
-          collocation_points_in_this_dim[index()[d]];
+    switch (mesh.basis(d)) {
+      case Spectral::Basis::SphericalHarmonic: {
+        switch (mesh.quadrature(d)) {
+          case Spectral::Quadrature::Gauss: {
+            const size_t n_theta = mesh.extents(d);
+            std::vector<double> theta(n_theta + 1);
+            DataVector temp(2 * n_theta + 1);
+            auto work = gsl::make_span(temp.data(), n_theta);
+            auto unused_weights =
+                gsl::make_span(temp.data() + n_theta, n_theta + 1);
+
+            int err = 0;
+            gaqd_(static_cast<int>(n_theta), theta.data(),
+                  unused_weights.data(), work.data(),
+                  static_cast<int>(unused_weights.size()), &err);
+            if (UNLIKELY(err != 0)) {
+              ERROR("gaqd error " << err << " in LogicalCoordinates");
+            }
+            for (IndexIterator<VolumeDim> index(mesh.extents()); index;
+                 ++index) {
+              logical_coords->get(d)[index.collapsed_index()] =
+                  theta[index()[d]];
+            }
+            break;
+          }
+          case Spectral::Quadrature::Equiangular: {
+            const size_t n_phi = mesh.extents(d);
+            const double two_pi_over_n_phi = 2.0 * M_PI / n_phi;
+            for (IndexIterator<VolumeDim> index(mesh.extents()); index;
+                 ++index) {
+              logical_coords->get(d)[index.collapsed_index()] =
+                  two_pi_over_n_phi * index()[d];
+            }
+            break;
+          }
+          default:
+            ERROR(
+                "Quadrature must be Gauss or Equiangular for Basis "
+                "SphericalHarmonic");
+        }
+        break;
+      }
+      // NOLINTNEXTLINE(bugprone-branch-clone)
+      case Spectral::Basis::Chebyshev:
+        [[fallthrough]];
+      case Spectral::Basis::Legendre:
+        [[fallthrough]];
+      case Spectral::Basis::FiniteDifference: {
+        const auto& collocation_points_in_this_dim =
+            Spectral::collocation_points(mesh.slice_through(d));
+        for (IndexIterator<VolumeDim> index(mesh.extents()); index; ++index) {
+          logical_coords->get(d)[index.collapsed_index()] =
+              collocation_points_in_this_dim[index()[d]];
+        }
+        break;
+      }
+      default:
+        ERROR("Missing basis case for logical_coordinates");
     }
   }
 }

src/NumericalAlgorithms/Spectral/LogicalCoordinates.hpp

@@ -37,10 +37,27 @@ struct not_null;
 /// @{
 /*!
  * \ingroup ComputationalDomainGroup
- * \brief Compute the logical coordinates in an Element.
+ * \brief Compute the logical coordinates of a Mesh in an Element.
  *
- * \details The logical coordinates are the collocation points associated to the
- * spectral basis functions and quadrature of the \p mesh.
+ * \details The logical coordinates are the collocation points
+ * associated with the Spectral::Basis and Spectral::Quadrature of the \p mesh.
+ * The Spectral::Basis determines the domain of the logical coordinates, while
+ * the Spectral::Quadrature determines their distribution.  For Legendre or
+ * Chebyshev bases, the logical coordinates are in the interval \f$[-1, 1]\f$.
+ * These bases may have either GaussLobatto or Gauss quadrature, which are not
+ * uniformly distributed, and either include (GaussLobatto) or do not include
+ * (Gauss) the end points.  For the FiniteDifference basis, the logical
+ * coordinates are again in the interval \f$[-1, 1]\f$.  This basis may have
+ * either FaceCentered or CellCentered quadrature, which are uniformly
+ * distributed, and either include (FaceCentered) or do not include
+ * (CellCentered) the end points.  The SphericalHarmonic basis corresponds to
+ * the spherical coordinates \f$(\theta, \phi)\f$ where the polar angle
+ * \f$\theta\f$ is in the interval \f$[0, \pi]\f$ and the azimuth \f$\phi\f$ is
+ * in the interval \f$[0, 2 \pi]\f$.  The polar angle has Gauss quadrature
+ * corresponding to the Legendre Gauss points of \f$\cos \theta\f$ (and thus
+ * have no points at the poles), while the azimuth has Equiangular quadrature
+ * which are distributed uniformly including the left endpoint, but not the
+ * right.
  *
  * \example
  * \snippet Test_LogicalCoordinates.cpp logical_coordinates_example

tests/Unit/NumericalAlgorithms/Spectral/CMakeLists.txt

@@ -32,5 +32,6 @@ target_link_libraries(
   ErrorHandling
   LinearOperators
   Spectral
+  SphericalHarmonics
   Utilities
   )

tests/Unit/NumericalAlgorithms/Spectral/Test_LogicalCoordinates.cpp

@@ -22,14 +22,37 @@
 #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp"
 #include "NumericalAlgorithms/Spectral/Mesh.hpp"
 #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
+#include "NumericalAlgorithms/SphericalHarmonics/Spherepack.hpp"
+
+namespace {
+
+void test_spherical_logical_coords() {
+  for (size_t l = 2; l < 5; ++l) {
+    const size_t nth = l + 1;
+    const size_t nph = 2 * l + 1;
+    const Mesh<2> mesh_s2{
+        {nth, nph},
+        {Spectral::Basis::SphericalHarmonic,
+         Spectral::Basis::SphericalHarmonic},
+        {Spectral::Quadrature::Gauss, Spectral::Quadrature::Equiangular}};
+    const auto xi = logical_coordinates(mesh_s2);
+    const ylm::Spherepack ylm(l, l);
+    const auto xi_expected = ylm.theta_phi_points();
+    CHECK(get<0>(xi) == xi_expected[0]);
+    CHECK(get<1>(xi) == xi_expected[1]);
+  }
+}
 
-namespace domain {
 SPECTRE_TEST_CASE("Unit.NumericalAlgorithms.Spectral.LogicalCoordinates",
                   "[Domain][Unit]") {
-  using Affine2d = CoordinateMaps::ProductOf2Maps<CoordinateMaps::Affine,
-                                                  CoordinateMaps::Affine>;
-  using Affine3d = CoordinateMaps::ProductOf3Maps<
-      CoordinateMaps::Affine, CoordinateMaps::Affine, CoordinateMaps::Affine>;
+  test_spherical_logical_coords();
+  using Affine2d =
+      domain::CoordinateMaps::ProductOf2Maps<domain::CoordinateMaps::Affine,
+                                             domain::CoordinateMaps::Affine>;
+  using Affine3d =
+      domain::CoordinateMaps::ProductOf3Maps<domain::CoordinateMaps::Affine,
+                                             domain::CoordinateMaps::Affine,
+                                             domain::CoordinateMaps::Affine>;
 
   const Mesh<1> mesh_1d{3, Spectral::Basis::Legendre,
                         Spectral::Quadrature::GaussLobatto};
@@ -41,20 +64,23 @@ SPECTRE_TEST_CASE("Unit.NumericalAlgorithms.Spectral.LogicalCoordinates",
                         Spectral::Basis::Legendre,
                         Spectral::Quadrature::GaussLobatto};
 
-  const CoordinateMaps::Affine x_map{-1.0, 1.0, -3.0, 7.0};
-  const CoordinateMaps::Affine y_map{-1.0, 1.0, -13.0, 47.0};
-  const CoordinateMaps::Affine z_map{-1.0, 1.0, -32.0, 74.0};
+  const domain::CoordinateMaps::Affine x_map{-1.0, 1.0, -3.0, 7.0};
+  const domain::CoordinateMaps::Affine y_map{-1.0, 1.0, -13.0, 47.0};
+  const domain::CoordinateMaps::Affine z_map{-1.0, 1.0, -32.0, 74.0};
 
-  const auto map_3d = make_coordinate_map<Frame::ElementLogical, Frame::Grid>(
-      Affine3d{x_map, y_map, z_map});
+  const auto map_3d =
+      domain::make_coordinate_map<Frame::ElementLogical, Frame::Grid>(
+          Affine3d{x_map, y_map, z_map});
 
   const auto x_3d = map_3d(logical_coordinates(mesh_3d));
   // [logical_coordinates_example]
 
-  const auto map_1d = make_coordinate_map<Frame::ElementLogical, Frame::Grid>(
-      CoordinateMaps::Affine{x_map});
-  const auto map_2d = make_coordinate_map<Frame::ElementLogical, Frame::Grid>(
-      Affine2d{x_map, y_map});
+  const auto map_1d =
+      domain::make_coordinate_map<Frame::ElementLogical, Frame::Grid>(
+          domain::CoordinateMaps::Affine{x_map});
+  const auto map_2d =
+      domain::make_coordinate_map<Frame::ElementLogical, Frame::Grid>(
+          Affine2d{x_map, y_map});
   const auto x_1d = map_1d(logical_coordinates(mesh_1d));
   const auto x_2d = map_2d(logical_coordinates(mesh_2d));
 
@@ -80,11 +106,11 @@ SPECTRE_TEST_CASE("Unit.NumericalAlgorithms.Spectral.LogicalCoordinates",
   CHECK(x_3d[2][0] == -32.0);
   CHECK(x_3d[2][15] == 74.0);
 
-  TestHelpers::db::test_compute_tag<Tags::LogicalCoordinates<1>>(
+  TestHelpers::db::test_compute_tag<domain::Tags::LogicalCoordinates<1>>(
       "ElementLogicalCoordinates");
-  TestHelpers::db::test_compute_tag<Tags::LogicalCoordinates<2>>(
+  TestHelpers::db::test_compute_tag<domain::Tags::LogicalCoordinates<2>>(
       "ElementLogicalCoordinates");
-  TestHelpers::db::test_compute_tag<Tags::LogicalCoordinates<3>>(
+  TestHelpers::db::test_compute_tag<domain::Tags::LogicalCoordinates<3>>(
       "ElementLogicalCoordinates");
 }
-}  // namespace domain
+}  // namespace