SwshCollocation.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstdlib>
7 #include <memory>
8 #include <sharp_cxx.h>
9 
10 #include "NumericalAlgorithms/Spectral/ComplexDataView.hpp"
12 #include "Utilities/Gsl.hpp"
13 
14 namespace Spectral {
15 namespace Swsh {
16 
17 /// \ingroup SwshGroup
18 /// \brief Convenience function for determining the number of spin-weighted
19 /// spherical harmonic collocation values that are stored for a given `l_max`
20 /// for a libsharp-compatible set of collocation points.
21 constexpr SPECTRE_ALWAYS_INLINE size_t
22 number_of_swsh_collocation_points(const size_t l_max) noexcept {
23  return (l_max + 1) * (2 * l_max + 1);
24 }
25 
26 // In the static caching mechanism, we permit an l_max up to this macro
27 // value. Higher l_max values may still be created manually using the
28 // `Collocation` constructor. If l_max's are used several times at higher value,
29 // consider increasing this value, but only after memory costs have been
30 // evaluated.
31 constexpr size_t collocation_maximum_l_max = 200;
32 
33 namespace detail {
34 // Helping functor to appropriately delete a stored `sharp_geom_info**`
35 struct DestroySharpGeometry {
36  void operator()(sharp_geom_info* to_delete) noexcept {
37  sharp_destroy_geom_info(to_delete);
38  }
39 };
40 } // namespace detail
41 
42 /// A container for reporting a single collocation point for libsharp compatible
43 /// data structures
45  size_t offset;
46  double theta;
47  double phi;
48 };
49 
50 /// \brief A wrapper class for the spherical harmonic library collocation data
51 ///
52 /// \details The currently chosen library for spin-weighted spherical harmonic
53 /// transforms is libsharp. The `Collocation` class stores the
54 /// libsharp `sharp_geom_info` object, which contains data about
55 /// 1. The angular collocation points used in spin-weighted spherical harmonic
56 /// transforms
57 /// 2. The memory representation of double-type values at those collocation
58 /// points
59 ///
60 /// \tparam Representation the ComplexRepresentation, either
61 /// `ComplexRepresentation::Interleaved` or
62 /// `ComplexRepresentation::RealsThenImags` compatible with the generated
63 /// `Collocation` - this is necessary because the stored libsharp type contains
64 /// memory stride information.
65 template <ComplexRepresentation Representation>
66 class Collocation {
67  public:
69  public:
70  /// create a new iterator. defaults to the start of the supplied object
72  const gsl::not_null<const Collocation<Representation>*> collocation,
73  const size_t start_index = 0) noexcept
74  : index_{start_index}, collocation_{collocation} {}
75 
76  /// recovers the data at the collocation point using a
77  /// `LibsharpCollocationPoint`, which stores the vector `offset` of the
78  /// location of the `theta`, `phi` point in libsharp compatible data
80  return LibsharpCollocationPoint{index_, collocation_->theta(index_),
81  collocation_->phi(index_)};
82  }
83  /// advance the iterator by one position (prefix)
85  ++index_;
86  return *this;
87  };
88  /// advance the iterator by one position (postfix)
89  // clang-tidy wants this to return a const iterator
90  CollocationConstIterator operator++(int)noexcept { // NOLINT
91  return CollocationConstIterator(collocation_, index_++);
92  }
93 
94  /// retreat the iterator by one position (prefix)
96  --index_;
97  return *this;
98  };
99  /// retreat the iterator by one position (prefix)
100  // clang-tidy wants this to return a const iterator
101  CollocationConstIterator operator--(int)noexcept { // NOLINT
102  return CollocationConstIterator(collocation_, index_--);
103  }
104 
105  // @{
106  /// (In)Equivalence checks both the object and index for the iterator
107  bool operator==(const CollocationConstIterator& rhs) const noexcept {
108  return index_ == rhs.index_ and collocation_ == rhs.collocation_;
109  }
110  bool operator!=(const CollocationConstIterator& rhs) const noexcept {
111  return not(*this == rhs);
112  }
113  // @}
114 
115  private:
116  size_t index_;
117  const Collocation<Representation>* const collocation_;
118  };
119 
120  /// The representation of the block of complex values, which sets the stride
121  /// inside the libsharp type
122  static constexpr ComplexRepresentation complex_representation =
123  Representation;
124 
125  /// \brief Generates the libsharp collocation information and stores it in
126  /// `geom_info_`.
127  ///
128  /// \note If you will potentially use the same `l_max` collocation set more
129  /// than once, it is probably better to use the
130  /// `precomputed_spherical_harmonic_collocation` function
131  explicit Collocation(size_t l_max) noexcept;
132 
133  /// default constructor required for iterator use
134  ~Collocation() = default;
135  Collocation() = default;
136  Collocation(const Collocation&) = default;
137  Collocation(Collocation&&) = default;
138  Collocation& operator=(Collocation&) = default;
139  Collocation& operator=(Collocation&&) = default;
140 
141  /// retrieve the `sharp_geom_info*` stored. This should largely be used only
142  /// for passing to other libsharp functions. Otherwise, access elements
143  /// through iterator or access functions.
144  sharp_geom_info* get_sharp_geom_info() const noexcept {
145  return geom_info_.get();
146  }
147 
148  /// Retrieve the \f$\theta\f$ value for a given index in a libsharp-compatible
149  /// array
150  double theta(size_t offset) const noexcept;
151  /// Retrieve the \f$\phi\f$ value for a given index in a libsharp-compatible
152  /// array
153  double phi(size_t offset) const noexcept;
154 
155  constexpr size_t l_max() const noexcept { return l_max_; }
156 
157  /// Compute the number of entries the libsharp-compatible data structure
158  /// should have
159  constexpr size_t size() const noexcept {
160  return (l_max_ + 1) * (2 * l_max_ + 1);
161  }
162 
163  // @{
164  /// Get a bidirectional iterator to the start of the grid. `operator*` for
165  /// that iterator gives a `LibsharpCollocationPoint` with members `offset`,
166  /// `theta`, and `phi` with operator *
168  return CollocationConstIterator{make_not_null(this), 0};
169  }
171  noexcept {
172  return begin();
173  }
174  // @}
175  // @{
176  /// Get a bidirectional iterator to the end of the grid. `operator*` for
177  /// that iterator gives a `LibsharpCollocationPoint` with members `offset`,
178  /// `theta`, and `phi` with operator *
180  return CollocationConstIterator{make_not_null(this), size()};
181  }
183  return end();
184  }
185  // @}
186 
187  private:
188  // an extra pointer layer is required for the peculiar way that libsharp
189  // constructs these values.
190  size_t l_max_ = 0;
192 };
193 
194 /// \brief precomputation function for those collocation grids that are
195 /// requested
196 ///
197 /// \details keeps a compile-time structure which acts as a thread-safe lookup
198 /// table for all l_max values that have been requested so far during execution,
199 /// so that the libsharp generation need not be re-run. If it has been
200 /// generated, it's returned by reference. Otherwise, the new grid is generated
201 /// and put in the lookup table before it is returned by reference.
202 template <ComplexRepresentation Representation>
204  size_t l_max) noexcept;
205 } // namespace Swsh
206 } // namespace Spectral
constexpr size_t size() const noexcept
Compute the number of entries the libsharp-compatible data structure should have. ...
Definition: SwshCollocation.hpp:159
constexpr size_t number_of_swsh_collocation_points(const size_t l_max) noexcept
Convenience function for determining the number of spin-weighted spherical harmonic collocation value...
Definition: SwshCollocation.hpp:22
LibsharpCollocationPoint operator*() const noexcept
recovers the data at the collocation point using a LibsharpCollocationPoint, which stores the vector ...
Definition: SwshCollocation.hpp:79
A wrapper class for the spherical harmonic library collocation data.
Definition: SwshCollocation.hpp:66
CollocationConstIterator operator++(int) noexcept
advance the iterator by one position (postfix)
Definition: SwshCollocation.hpp:90
sharp_geom_info * get_sharp_geom_info() const noexcept
retrieve the sharp_geom_info* stored. This should largely be used only for passing to other libsharp ...
Definition: SwshCollocation.hpp:144
Definition: Determinant.hpp:11
#define SPECTRE_ALWAYS_INLINE
Always inline a function. Only use this if you benchmarked the code.
Definition: ForceInline.hpp:20
bool operator!=(const CollocationConstIterator &rhs) const noexcept
(In)Equivalence checks both the object and index for the iterator
Definition: SwshCollocation.hpp:110
tnsr::iaa< DataType, SpatialDim, Frame > phi(const Scalar< DataType > &lapse, const tnsr::i< DataType, SpatialDim, Frame > &deriv_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::iJ< DataType, SpatialDim, Frame > &deriv_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ijj< DataType, SpatialDim, Frame > &deriv_spatial_metric) noexcept
Computes the auxiliary variable used by the generalized harmonic formulation of Einstein&#39;s equations...
Definition: ComputeGhQuantities.cpp:22
A container for reporting a single collocation point for libsharp compatible data structures...
Definition: SwshCollocation.hpp:44
CollocationConstIterator operator--(int) noexcept
retreat the iterator by one position (prefix)
Definition: SwshCollocation.hpp:101
CollocationConstIterator(const gsl::not_null< const Collocation< Representation > *> collocation, const size_t start_index=0) noexcept
create a new iterator. defaults to the start of the supplied object
Definition: SwshCollocation.hpp:71
Defines macro to always inline a function.
CollocationConstIterator & operator--() noexcept
retreat the iterator by one position (prefix)
Definition: SwshCollocation.hpp:95
bool operator==(const CollocationConstIterator &rhs) const noexcept
(In)Equivalence checks both the object and index for the iterator
Definition: SwshCollocation.hpp:107
Collocation< Representation >::CollocationConstIterator cbegin() const noexcept
Get a bidirectional iterator to the start of the grid. operator* for that iterator gives a LibsharpCo...
Definition: SwshCollocation.hpp:170
Defines functions and classes from the GSL.
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion, but it may be necessary to perform the conversion explicitly when type deduction is desired.
Definition: Gsl.hpp:863
Functionality associated with a particular choice of basis functions and quadrature for spectral oper...
Definition: Chebyshev.cpp:19
ComplexRepresentation
A set of labels for the possible representations of complex numbers for the ComplexDataView and the c...
Definition: ComplexDataView.hpp:57
Collocation< Representation >::CollocationConstIterator begin() const noexcept
Get a bidirectional iterator to the start of the grid. operator* for that iterator gives a LibsharpCo...
Definition: SwshCollocation.hpp:167
Collocation< Representation >::CollocationConstIterator end() const noexcept
Get a bidirectional iterator to the end of the grid. operator* for that iterator gives a LibsharpColl...
Definition: SwshCollocation.hpp:179
Collocation< Representation >::CollocationConstIterator cend() const noexcept
Get a bidirectional iterator to the end of the grid. operator* for that iterator gives a LibsharpColl...
Definition: SwshCollocation.hpp:182
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12
CollocationConstIterator & operator++() noexcept
advance the iterator by one position (prefix)
Definition: SwshCollocation.hpp:84
const Collocation< Representation > & precomputed_collocation(const size_t l_max) noexcept
precomputation function for those collocation grids that are requested
Definition: SwshCollocation.cpp:97