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/ComplexDataView.hpp"
      13             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      14             : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
      15             : #include "NumericalAlgorithms/Spectral/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