Line data Source code
1 0 : // Distributed under the MIT License. 2 : // See LICENSE.txt for details. 3 : 4 : #pragma once 5 : 6 : #include <memory> 7 : #include <sharp_cxx.h> 8 : 9 : #include "DataStructures/DataVector.hpp" 10 : #include "DataStructures/Tensor/Tensor.hpp" 11 : #include "NumericalAlgorithms/Spectral/Basis.hpp" 12 : #include "NumericalAlgorithms/Spectral/Mesh.hpp" 13 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp" 14 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/ComplexDataView.hpp" 15 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshSettings.hpp" 16 : #include "Utilities/ForceInline.hpp" 17 : #include "Utilities/Gsl.hpp" 18 : 19 : namespace Spectral { 20 : namespace Swsh { 21 : 22 : /// \ingroup SwshGroup 23 : /// \brief Convenience function for determining the number of spin-weighted 24 : /// spherical harmonic collocation values that are stored for a given `l_max` 25 : /// for a libsharp-compatible set of collocation points. 26 : constexpr SPECTRE_ALWAYS_INLINE size_t 27 1 : number_of_swsh_collocation_points(const size_t l_max) { 28 : return (l_max + 1) * (2 * l_max + 1); 29 : } 30 : 31 : /// \ingroup SwshGroup 32 : /// \brief Returns the number of spin-weighted spherical harmonic collocation 33 : /// values in \f$\theta\f$ for a libsharp-compatible set of collocation 34 : /// points. 35 : /// 36 : /// \details The full number of collocation points is the product of the number 37 : /// of \f$\theta\f$ points and the number of \f$\phi\f$ points (a 'rectangular' 38 : /// grid). 39 : constexpr SPECTRE_ALWAYS_INLINE size_t 40 1 : number_of_swsh_theta_collocation_points(const size_t l_max) { 41 : return (l_max + 1); 42 : } 43 : 44 : /// \ingroup SwshGroup 45 : /// \brief Returns the number of spin-weighted spherical harmonic collocation 46 : /// values in \f$\phi\f$ for a libsharp-compatible set of collocation 47 : /// points. 48 : /// 49 : /// \details The full number of collocation points is the product of the number 50 : /// of \f$\theta\f$ points and the number of \f$\phi\f$ points (a 'rectangular' 51 : /// grid). 52 : constexpr SPECTRE_ALWAYS_INLINE size_t 53 1 : number_of_swsh_phi_collocation_points(const size_t l_max) { 54 : return (2 * l_max + 1); 55 : } 56 : 57 : /// \ingroup SwshGroup 58 : /// \brief Obtain the three-dimensional mesh associated with a 59 : /// libsharp-compatible sequence of spherical nodal shells. 60 : /// \warning This is to be used only for operations in the radial direction! The 61 : /// angular collocation points should not be regarded as having a useful product 62 : /// representation for e.g. taking derivatives in just the theta direction. 63 : /// Instead, use spin-weighted utilities for all angular operations. 64 1 : SPECTRE_ALWAYS_INLINE Mesh<3> swsh_volume_mesh_for_radial_operations( 65 : const size_t l_max, const size_t number_of_radial_points) { 66 : return Mesh<3>{{{number_of_swsh_phi_collocation_points(l_max), 67 : number_of_swsh_theta_collocation_points(l_max), 68 : number_of_radial_points}}, 69 : Spectral::Basis::Legendre, 70 : Spectral::Quadrature::GaussLobatto}; 71 : } 72 : 73 : namespace detail { 74 : // Helping functor to appropriately delete a stored `sharp_geom_info**` 75 : struct DestroySharpGeometry { 76 : void operator()(sharp_geom_info* to_delete) { 77 : sharp_destroy_geom_info(to_delete); 78 : } 79 : }; 80 : } // namespace detail 81 : 82 : /// A container for reporting a single collocation point for libsharp compatible 83 : /// data structures 84 1 : struct LibsharpCollocationPoint { 85 0 : size_t offset; 86 0 : double theta; 87 0 : double phi; 88 : }; 89 : 90 : /// \brief A wrapper class for the spherical harmonic library collocation data 91 : /// 92 : /// \details The currently chosen library for spin-weighted spherical harmonic 93 : /// transforms is libsharp. The `CollocationMetadata` class stores the 94 : /// libsharp `sharp_geom_info` object, which contains data about 95 : /// 1. The angular collocation points used in spin-weighted spherical harmonic 96 : /// transforms 97 : /// 2. The memory representation of double-type values at those collocation 98 : /// points 99 : /// 100 : /// \tparam Representation the ComplexRepresentation, either 101 : /// `ComplexRepresentation::Interleaved` or 102 : /// `ComplexRepresentation::RealsThenImags` compatible with the generated 103 : /// `CollocationMetadata` - this is necessary because the stored libsharp type 104 : /// contains memory stride information. 105 : template <ComplexRepresentation Representation> 106 1 : class CollocationMetadata { 107 : public: 108 0 : class CollocationConstIterator { 109 : public: 110 : /// create a new iterator. defaults to the start of the supplied object 111 1 : explicit CollocationConstIterator( 112 : const gsl::not_null<const CollocationMetadata<Representation>*> 113 : collocation, 114 : const size_t start_index = 0) 115 : : index_{start_index}, collocation_{collocation} {} 116 : 117 : /// recovers the data at the collocation point using a 118 : /// `LibsharpCollocationPoint`, which stores the vector `offset` of the 119 : /// location of the `theta`, `phi` point in libsharp compatible data 120 1 : LibsharpCollocationPoint operator*() const { 121 : return LibsharpCollocationPoint{index_, collocation_->theta(index_), 122 : collocation_->phi(index_)}; 123 : } 124 : /// advance the iterator by one position (prefix) 125 1 : CollocationConstIterator& operator++() { 126 : ++index_; 127 : return *this; 128 : }; 129 : /// advance the iterator by one position (postfix) 130 : // clang-tidy wants this to return a const iterator 131 1 : CollocationConstIterator operator++(int) { // NOLINT 132 : return CollocationConstIterator(collocation_, index_++); 133 : } 134 : 135 : /// retreat the iterator by one position (prefix) 136 1 : CollocationConstIterator& operator--() { 137 : --index_; 138 : return *this; 139 : }; 140 : /// retreat the iterator by one position (prefix) 141 : // clang-tidy wants this to return a const iterator 142 1 : CollocationConstIterator operator--(int) { // NOLINT 143 : return CollocationConstIterator(collocation_, index_--); 144 : } 145 : 146 : /// @{ 147 : /// (In)Equivalence checks both the object and index for the iterator 148 1 : bool operator==(const CollocationConstIterator& rhs) const { 149 : return index_ == rhs.index_ and collocation_ == rhs.collocation_; 150 : } 151 1 : bool operator!=(const CollocationConstIterator& rhs) const { 152 : return not(*this == rhs); 153 : } 154 : /// @} 155 : 156 : private: 157 0 : size_t index_; 158 0 : const CollocationMetadata<Representation>* const collocation_; 159 : }; 160 : 161 : /// The representation of the block of complex values, which sets the stride 162 : /// inside the libsharp type 163 1 : static constexpr ComplexRepresentation complex_representation = 164 : Representation; 165 : 166 : /// \brief Generates the libsharp collocation information and stores it in 167 : /// `geom_info_`. 168 : /// 169 : /// \note If you will potentially use the same `l_max` collocation set more 170 : /// than once, it is probably better to use the 171 : /// `precomputed_spherical_harmonic_collocation` function 172 1 : explicit CollocationMetadata(size_t l_max); 173 : 174 : /// default constructor required for iterator use 175 1 : ~CollocationMetadata() = default; 176 0 : CollocationMetadata() = default; 177 0 : CollocationMetadata(const CollocationMetadata&) = default; 178 0 : CollocationMetadata(CollocationMetadata&&) = default; 179 0 : CollocationMetadata& operator=(CollocationMetadata&) = default; 180 0 : CollocationMetadata& operator=(CollocationMetadata&&) = default; 181 : 182 : /// retrieve the `sharp_geom_info*` stored. This should largely be used only 183 : /// for passing to other libsharp functions. Otherwise, access elements 184 : /// through iterator or access functions. 185 1 : sharp_geom_info* get_sharp_geom_info() const { return geom_info_.get(); } 186 : 187 : /// Retrieve the \f$\theta\f$ value for a given index in a libsharp-compatible 188 : /// array 189 1 : double theta(size_t offset) const; 190 : /// Retrieve the \f$\phi\f$ value for a given index in a libsharp-compatible 191 : /// array 192 1 : double phi(size_t offset) const; 193 : 194 0 : constexpr size_t l_max() const { return l_max_; } 195 : 196 : /// Compute the number of entries the libsharp-compatible data structure 197 : /// should have 198 1 : constexpr size_t size() const { return (l_max_ + 1) * (2 * l_max_ + 1); } 199 : 200 : /// @{ 201 : /// Get a bidirectional iterator to the start of the grid. `operator*` for 202 : /// that iterator gives a `LibsharpCollocationPoint` with members `offset`, 203 : /// `theta`, and `phi` 204 1 : CollocationMetadata<Representation>::CollocationConstIterator begin() const { 205 : return CollocationConstIterator{make_not_null(this), 0}; 206 : } 207 1 : CollocationMetadata<Representation>::CollocationConstIterator cbegin() const { 208 : return begin(); 209 : } 210 : /// @} 211 : /// @{ 212 : /// Get a bidirectional iterator to the end of the grid. `operator*` for 213 : /// that iterator gives a `LibsharpCollocationPoint` with members `offset`, 214 : /// `theta`, and `phi` 215 1 : CollocationMetadata<Representation>::CollocationConstIterator end() const { 216 : return CollocationConstIterator{make_not_null(this), size()}; 217 : } 218 1 : CollocationMetadata<Representation>::CollocationConstIterator cend() const { 219 : return end(); 220 : } 221 : /// @} 222 : 223 : private: 224 : // an extra pointer layer is required for the peculiar way that libsharp 225 : // constructs these values. 226 0 : size_t l_max_ = 0; 227 0 : std::unique_ptr<sharp_geom_info, detail::DestroySharpGeometry> geom_info_; 228 : }; 229 : 230 : /// \brief precomputation function for those collocation grids that are 231 : /// requested 232 : /// 233 : /// \details keeps a compile-time structure which acts as a thread-safe lookup 234 : /// table for all l_max values that have been requested so far during execution, 235 : /// so that the libsharp generation need not be re-run. If it has been 236 : /// generated, it's returned by reference. Otherwise, the new grid is generated 237 : /// and put in the lookup table before it is returned by reference. 238 : template <ComplexRepresentation Representation> 239 1 : const CollocationMetadata<Representation>& cached_collocation_metadata( 240 : size_t l_max); 241 : 242 : /// \brief Store the libsharp-compatible collocation grid and corresponding 243 : /// unit-sphere cartesian grid in the supplied buffers. 244 1 : void create_angular_and_cartesian_coordinates( 245 : const gsl::not_null<tnsr::i<DataVector, 3>*> cartesian_coordinates, 246 : const gsl::not_null< 247 : tnsr::i<DataVector, 2, ::Frame::Spherical<::Frame::Inertial>>*> 248 : angular_coordinates, 249 : size_t l_max); 250 : } // namespace Swsh 251 : } // namespace Spectral