From e4c25937a8db6fe631b1e68032e5b17d01c6cfbf Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Fri, 20 Oct 2023 22:56:02 -0700 Subject: [PATCH] Enable Ylm meshes --- .../Spectral/CMakeLists.txt | 1 + src/NumericalAlgorithms/Spectral/Mesh.cpp | 8 --- src/NumericalAlgorithms/Spectral/Spectral.cpp | 26 +++++++- .../Spectral/SphericalHarmonics.cpp | 60 +++++++++++++++++++ 4 files changed, 84 insertions(+), 11 deletions(-) create mode 100644 src/NumericalAlgorithms/Spectral/SphericalHarmonics.cpp diff --git a/src/NumericalAlgorithms/Spectral/CMakeLists.txt b/src/NumericalAlgorithms/Spectral/CMakeLists.txt index 07d9fe911acc6..5b5a2373191ac 100644 --- a/src/NumericalAlgorithms/Spectral/CMakeLists.txt +++ b/src/NumericalAlgorithms/Spectral/CMakeLists.txt @@ -18,6 +18,7 @@ spectre_target_sources( Projection.cpp Quadrature.cpp Spectral.cpp + SphericalHarmonics.cpp ) spectre_target_headers( diff --git a/src/NumericalAlgorithms/Spectral/Mesh.cpp b/src/NumericalAlgorithms/Spectral/Mesh.cpp index fb213b59ea9e1..0e80e29f3e2d8 100644 --- a/src/NumericalAlgorithms/Spectral/Mesh.cpp +++ b/src/NumericalAlgorithms/Spectral/Mesh.cpp @@ -33,8 +33,6 @@ template Mesh::Mesh(const size_t isotropic_extents, const Spectral::Basis basis, const Spectral::Quadrature quadrature) { static_assert(sizeof(Mesh) == (Dim == 0 ? 2 : 2 * Dim)); - ASSERT(basis != Spectral::Basis::SphericalHarmonic, - "SphericalHarmonic is not a valid basis for the Mesh"); ASSERT(isotropic_extents <= 255, "Cannot have more than 255 grid points"); extents_[0] = isotropic_extents; quadrature_and_basis_[0] = combine(basis, quadrature); @@ -52,8 +50,6 @@ template Mesh::Mesh(const std::array& extents, const Spectral::Basis basis, const Spectral::Quadrature quadrature) { - ASSERT(basis != Spectral::Basis::SphericalHarmonic, - "SphericalHarmonic is not a valid basis for the Mesh"); if constexpr (Dim > 0) { ASSERT( extents[0] <= 255, @@ -86,10 +82,6 @@ template Mesh::Mesh(const std::array& extents, const std::array& bases, const std::array& quadratures) { - for (auto it = bases.begin(); it != bases.end(); it++) { - ASSERT(*it != Spectral::Basis::SphericalHarmonic, - "SphericalHarmonic is not a valid basis for the Mesh"); - } if constexpr (Dim > 0) { ASSERT( extents[0] <= 255, diff --git a/src/NumericalAlgorithms/Spectral/Spectral.cpp b/src/NumericalAlgorithms/Spectral/Spectral.cpp index 2b0c446b08cdf..7e4de7772ac48 100644 --- a/src/NumericalAlgorithms/Spectral/Spectral.cpp +++ b/src/NumericalAlgorithms/Spectral/Spectral.cpp @@ -428,6 +428,8 @@ struct DifferentiationMatrixGenerator { diff_matrix(i, j) *= inv_delta; } } + } else if constexpr (BasisType == Basis::SphericalHarmonic) { + ERROR("Not implemented"); } else { const DataVector& bary_weights = barycentric_weights(num_points); @@ -797,9 +799,7 @@ template decltype(auto) get_spectral_quantity_for_mesh(F&& f, const Mesh<1>& mesh) { const auto num_points = mesh.extents(0); // Switch on runtime values of basis and quadrature to select - // corresponding template specialization. For basis functions spanning - // multiple dimensions we can generalize this function to take a - // higher-dimensional Mesh. + // corresponding template specialization. switch (mesh.basis(0)) { case Basis::Legendre: switch (mesh.quadrature(0)) { @@ -829,6 +829,20 @@ decltype(auto) get_spectral_quantity_for_mesh(F&& f, const Mesh<1>& mesh) { default: ERROR("Missing quadrature case for spectral quantity"); } + case Basis::SphericalHarmonic: + switch (mesh.quadrature(0)) { + case Quadrature::Gauss: // [0, pi] direction + return f(std::integral_constant{}, + std::integral_constant{}, + num_points); + case Quadrature::Equiangular: // [0, 2 pi) direction + return f( + std::integral_constant{}, + std::integral_constant{}, + num_points); + default: + ERROR("Missing quadrature case for spectral quantity"); + } case Basis::FiniteDifference: switch (mesh.quadrature(0)) { case Quadrature::CellCentered: @@ -932,6 +946,12 @@ GENERATE_INSTANTIATIONS(INSTANTIATE, #undef QUAD #undef INSTANTIATE +template const DataVector& Spectral::collocation_points< + Spectral::Basis::SphericalHarmonic, Spectral::Quadrature::Gauss>(size_t); +template const DataVector& + Spectral::collocation_points(size_t); + template const DataVector& Spectral::collocation_points(size_t); diff --git a/src/NumericalAlgorithms/Spectral/SphericalHarmonics.cpp b/src/NumericalAlgorithms/Spectral/SphericalHarmonics.cpp new file mode 100644 index 0000000000000..8463c94965e5d --- /dev/null +++ b/src/NumericalAlgorithms/Spectral/SphericalHarmonics.cpp @@ -0,0 +1,60 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#include "NumericalAlgorithms/Spectral/Spectral.hpp" + +#include + +#include "DataStructures/DataVector.hpp" +#include "DataStructures/Matrix.hpp" +#include "Utilities/ErrorHandling/Error.hpp" + +namespace Spectral { + +// Algorithms to compute spherical harmonic basis functions +// These functions specialize the templates declared in `Spectral.hpp`. + +template <> +std::pair compute_collocation_points_and_weights< + Basis::SphericalHarmonic, Quadrature::Gauss>(const size_t num_points) { + return compute_collocation_points_and_weights(num_points); +} + +template <> +std::pair compute_collocation_points_and_weights< + Basis::SphericalHarmonic, Quadrature::Equiangular>( + const size_t num_points) { + return compute_collocation_points_and_weights( + num_points); +} + +template <> +DataVector compute_basis_function_value( + const size_t /*k*/, const DataVector& /*x*/) { + ERROR("not implemented"); +} + +template <> +DataVector compute_inverse_weight_function_values( + const DataVector& /*x*/) { + ERROR("not implemented"); +} + +template <> +double compute_basis_function_normalization_square( + const size_t /*k*/) { + ERROR("not implemented"); +} + +template +Matrix spectral_indefinite_integral_matrix(size_t num_points); + +template <> +Matrix spectral_indefinite_integral_matrix( + const size_t /*num_points*/) { + ERROR("not implemented"); +} + +} // namespace Spectral