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
|
