Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details
3 :
4 : #pragma once
5 :
6 : #include <complex>
7 : #include <memory>
8 : #include <sharp_cxx.h>
9 :
10 : #include "DataStructures/ComplexModalVector.hpp"
11 : #include "DataStructures/SpinWeighted.hpp"
12 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshCollocation.hpp"
13 : #include "NumericalAlgorithms/SpinWeightedSphericalHarmonics/SwshSettings.hpp"
14 : #include "Utilities/ConstantExpressions.hpp"
15 : #include "Utilities/ForceInline.hpp"
16 : #include "Utilities/Gsl.hpp"
17 :
18 : namespace Spectral {
19 : namespace Swsh {
20 :
21 : /// \ingroup SwshGroup
22 : /// \brief Convenience function for determining the number of spin-weighted
23 : /// spherical harmonics coefficients that are stored for a given `l_max`
24 : ///
25 : /// \details This includes the factor of 2 associated with needing
26 : /// to store both the transform of the real and imaginary parts, so is the
27 : /// full size of the result of a libsharp swsh transform.
28 : ///
29 : /// \note Assumes the triangular libsharp representation is used.
30 : constexpr SPECTRE_ALWAYS_INLINE size_t
31 1 : size_of_libsharp_coefficient_vector(const size_t l_max) {
32 : return (l_max + 1) * (l_max + 2); // "triangular" representation
33 : }
34 :
35 : /*!
36 : * \ingroup SwshGroup
37 : * \brief Compute the relative sign change necessary to convert between the
38 : * libsharp basis for spin weight `from_spin_weight` to the basis for spin
39 : * weight `to_spin_weight`, for the real component coefficients if `real` is
40 : * true, otherwise for the imaginary component coefficients. The sign change for
41 : * a given coefficient is equivalent to the product of
42 : * `sharp_swsh_sign(from_spin, m, real) * sharp_swsh_sign(to_spin, m,
43 : * real)`. Due to the form of the signs, it does not end up depending on m (the
44 : * m's in the power of \f$-1\f$'s cancel).
45 : * For full details of the libsharp sign conventions, see the documentation for
46 : * TransformJob.
47 : *
48 : * \details The sign change is obtained by the
49 : * difference between the libsharp convention and the convention which uses:
50 : *
51 : * \f[
52 : * {}_s Y_{\ell m} (\theta, \phi) = (-1)^m \sqrt{\frac{2 l + 1}{4 \pi}}
53 : * D^{\ell}{}_{-m s}(\phi, \theta, 0).
54 : * \f]
55 : *
56 : * See \cite Goldberg1966uu
57 : */
58 1 : constexpr SPECTRE_ALWAYS_INLINE double sharp_swsh_sign_change(
59 : const int from_spin_weight, const int to_spin_weight, const bool real) {
60 : if (real) {
61 : return (from_spin_weight == 0 ? -1.0 : 1.0) *
62 : (from_spin_weight >= 0 ? -1.0 : 1.0) *
63 : ((from_spin_weight < 0 and (from_spin_weight % 2 == 0)) ? -1.0
64 : : 1.0) *
65 : (to_spin_weight == 0 ? -1.0 : 1.0) *
66 : (to_spin_weight >= 0 ? -1.0 : 1.0) *
67 : ((to_spin_weight < 0 and (to_spin_weight % 2 == 0)) ? -1.0 : 1.0);
68 : }
69 : return (from_spin_weight == 0 ? -1.0 : 1.0) *
70 : ((from_spin_weight < 0 and (from_spin_weight % 2 == 0)) ? 1.0 : -1.0) *
71 : (to_spin_weight == 0 ? -1.0 : 1.0) *
72 : ((to_spin_weight < 0 and (to_spin_weight % 2 == 0)) ? 1.0 : -1.0);
73 : }
74 :
75 : /*!
76 : * \ingroup SwshGroup
77 : * \brief Compute the sign change between the libsharp convention and the set
78 : * of spin-weighted spherical harmonics given by the relation to the Wigner
79 : * rotation matrices.
80 : *
81 : * \details The sign change is obtained via the difference between the libsharp
82 : * convention and the convention which uses:
83 : *
84 : * \f[
85 : * {}_s Y_{\ell m} (\theta, \phi) = (-1)^m \sqrt{\frac{2 l + 1}{4 \pi}}
86 : * D^{\ell}{}_{-m s}(\phi, \theta, 0).
87 : * \f]
88 : *
89 : * See \cite Goldberg1966uu.
90 : * The sign change is computed for the real component coefficients if `real` is
91 : * true, otherwise for the imaginary component coefficients. For full details on
92 : * the sign convention used in libsharp, see the documentation for TransformJob.
93 : * This function outputs the \f$\mathrm{sign}(m, s, \mathrm{real})\f$ necessary
94 : * to produce the conversion between Goldberg moments and libsharp moments:
95 : *
96 : * \f{align*}{
97 : * {}_s Y_{\ell m}^{\mathrm{real}, \mathrm{sharp}} =& \mathrm{sign}(m, s,
98 : * \mathrm{real=true}) {}_s Y_{\ell m}^{\mathrm{Goldberg}}\\
99 : * {}_s Y_{\ell m}^{\mathrm{imag}, \mathrm{sharp}} =&
100 : * \mathrm{sign}(m, s, \mathrm{real=false}) {}_s Y_{\ell
101 : * m}^{\mathrm{Goldberg}}.
102 : * \f}
103 : *
104 : * Note that in this equation, the "real" and "imag" superscripts refer to the
105 : * set of basis functions used for the decomposition of the real and imaginary
106 : * part of the spin-weighted collocation points, not real or imaginary parts of
107 : * the basis functions themselves.
108 : */
109 1 : constexpr SPECTRE_ALWAYS_INLINE double sharp_swsh_sign(const int spin_weight,
110 : const int m,
111 : const bool real) {
112 : if (real) {
113 : if (m >= 0) {
114 : return (spin_weight > 0 ? -1.0 : 1.0) *
115 : ((spin_weight < 0 and (spin_weight % 2 == 0)) ? -1.0 : 1.0);
116 : }
117 : return (spin_weight == 0 ? -1.0 : 1.0) *
118 : ((spin_weight >= 0 and (m % 2 == 0)) ? -1.0 : 1.0) *
119 : (spin_weight < 0 and (m + spin_weight) % 2 == 0 ? -1.0 : 1.0);
120 : }
121 : if (m >= 0) {
122 : return (spin_weight == 0 ? -1.0 : 1.0) *
123 : ((spin_weight < 0 and (spin_weight % 2 == 0)) ? 1.0 : -1.0);
124 : }
125 : return (spin_weight == 0 ? -1.0 : 1.0) *
126 : ((spin_weight >= 0 and (m % 2 == 0)) ? -1.0 : 1.0) *
127 : ((spin_weight < 0 and (m + spin_weight) % 2 != 0) ? -1.0 : 1.0);
128 : }
129 :
130 : namespace detail {
131 : struct DestroySharpAlm {
132 : void operator()(sharp_alm_info* to_delete) {
133 : sharp_destroy_alm_info(to_delete);
134 : }
135 : };
136 : } // namespace detail
137 :
138 : /// Points to a single pair of modes in a libsharp-compatible
139 : /// spin-weighted spherical harmonic modal representation.
140 1 : struct LibsharpCoefficientInfo {
141 0 : size_t transform_of_real_part_offset;
142 0 : size_t transform_of_imag_part_offset;
143 0 : size_t l_max;
144 0 : size_t l;
145 0 : size_t m;
146 : };
147 :
148 : /*!
149 : * \ingroup SwshGroup
150 : * \brief A container for libsharp metadata for the spin-weighted spherical
151 : * harmonics modal representation.
152 : *
153 : * \details
154 : * The CoefficientsMetadata class acts as a memory-safe container for a
155 : * `sharp_alm_info*`, required for use of libsharp transform utilities.
156 : * The libsharp utilities are currently constructed to only provide user
157 : * functions with collocation data for spin-weighted functions and
158 : * derivatives. This class also provides an iterator for
159 : * easily traversing a libsharp-compatible modal representation.
160 : *
161 : * \note The libsharp representation of coefficients is altered from the
162 : * standard mathematical definitions in a nontrivial way. There are a number of
163 : * important features to the data storage of the coefficients.
164 : * - they are stored as a set of complex values, but each vector of complex
165 : * values is the transform of only the real or imaginary part of the
166 : * collocation data
167 : * - because each vector of complex coefficients is related to the transform of
168 : * a set of doubles, only (about) half of the m's are stored (m >= 0), because
169 : * the remaining m modes are determinable by conjugation from the positive m
170 : * modes, given that they represent the transform of a purely real or purely
171 : * imaginary collocation quantity
172 : * - they are stored in an l-varies-fastest, triangular representation. To be
173 : * concrete, for an l_max=2, the order of coefficient storage is (l, m):
174 : * [(0, 0), (1, 0), (2, 0), (1, 1), (2, 1), (2, 2)]
175 : * - due to the restriction of representing only the transform of real
176 : * quantities, the m=0 modes always have vanishing imaginary component.
177 : */
178 1 : class CoefficientsMetadata {
179 : public:
180 : /// An iterator for easily traversing a libsharp-compatible spin-weighted
181 : /// spherical harmonic modal representation.
182 : /// The `operator*()` returns a `LibsharpCoefficientInfo`, which contains two
183 : /// offsets, `transform_of_real_part_offset` and
184 : /// `transform_of_imag_part_offset`, and the `l_max`, `l` and `m` associated
185 : /// with the values at those offsets.
186 : ///
187 : /// \note this currently assumes, as do many of the utilities in this file,
188 : /// that the libsharp representation is chosen to be the triangular
189 : /// coefficient representation. If alternative representations are desired,
190 : /// alterations will be needed.
191 1 : class CoefficientsIndexIterator {
192 : public:
193 0 : explicit CoefficientsIndexIterator(const size_t l_max,
194 : const size_t start_l = 0,
195 : const size_t start_m = 0)
196 : : m_{start_m}, l_{start_l}, l_max_{l_max} {}
197 :
198 0 : LibsharpCoefficientInfo operator*() const {
199 : // permit dereferencing the iterator only if the return represents a
200 : // viable location in the coefficient vector
201 : ASSERT(l_ <= l_max_ && m_ <= l_max_, "coefficients iterator overflow");
202 : size_t offset = ((3 + 2 * l_max_ - m_) * m_) / 2 + l_ - m_;
203 : return LibsharpCoefficientInfo{
204 : offset,
205 : offset +
206 : Spectral::Swsh::size_of_libsharp_coefficient_vector(l_max_) / 2,
207 : l_max_, l_, m_};
208 : }
209 :
210 : /// advance the iterator by one position (prefix)
211 1 : CoefficientsIndexIterator& operator++() {
212 : // permit altering the iterator only if the new value is between the
213 : // anticipated begin and end, inclusive
214 : ASSERT(l_ <= l_max_ && m_ <= l_max_, "coefficients iterator overflow");
215 : if (l_ == l_max_) {
216 : ++m_;
217 : l_ = m_;
218 : } else {
219 : ++l_;
220 : }
221 : return *this;
222 : };
223 :
224 : /// advance the iterator by one position (postfix)
225 : // NOLINTNEXTLINE(readability-const-return-type)
226 1 : const CoefficientsIndexIterator operator++(int) {
227 : auto pre_increment = *this;
228 : ++*this;
229 : return pre_increment;
230 : }
231 :
232 : /// retreat the iterator by one position (prefix)
233 1 : CoefficientsIndexIterator& operator--() {
234 : // permit altering the iterator only if the new value is between the
235 : // anticipated begin and end, inclusive
236 : ASSERT(l_ <= l_max_ + 1 && m_ <= l_max_ + 1,
237 : "coefficients iterator overflow");
238 : if (l_ == m_) {
239 : --m_;
240 : l_ = l_max_;
241 : } else {
242 : --l_;
243 : }
244 : return *this;
245 : }
246 :
247 : /// retreat the iterator by one position (postfix)
248 : // NOLINTNEXTLINE(readability-const-return-type)
249 1 : const CoefficientsIndexIterator operator--(int) {
250 : auto pre_decrement = *this;
251 : --*this;
252 : return pre_decrement;
253 : }
254 :
255 : /// @{
256 : /// (In)Equivalence checks the object as well as the l and m current
257 : /// position.
258 1 : bool operator==(const CoefficientsIndexIterator& rhs) const {
259 : return m_ == rhs.m_ and l_ == rhs.l_ and l_max_ == rhs.l_max_;
260 : }
261 1 : bool operator!=(const CoefficientsIndexIterator& rhs) const {
262 : return not(*this == rhs);
263 : }
264 : /// @}
265 :
266 : private:
267 0 : size_t m_;
268 0 : size_t l_;
269 0 : size_t l_max_;
270 : };
271 :
272 0 : explicit CoefficientsMetadata(size_t l_max);
273 :
274 0 : ~CoefficientsMetadata() = default;
275 0 : CoefficientsMetadata() = default;
276 0 : CoefficientsMetadata(const CoefficientsMetadata&) = delete;
277 0 : CoefficientsMetadata(CoefficientsMetadata&&) = default;
278 0 : CoefficientsMetadata& operator=(const CoefficientsMetadata&) = delete;
279 0 : CoefficientsMetadata& operator=(CoefficientsMetadata&&) = default;
280 0 : sharp_alm_info* get_sharp_alm_info() const { return alm_info_.get(); }
281 :
282 0 : size_t l_max() const { return l_max_; }
283 :
284 : /// returns the number of (complex) entries in a libsharp-compatible
285 : /// coefficients vector. This includes the factor of 2 associated with needing
286 : /// to store both the transform of the real and imaginary parts, so is the
287 : /// full size of the result of a libsharp swsh transform.
288 1 : size_t size() const { return size_of_libsharp_coefficient_vector(l_max_); }
289 :
290 : /// @{
291 : /// \brief Get a bidirectional iterator to the start of the series of modes.
292 1 : CoefficientsMetadata::CoefficientsIndexIterator begin() const {
293 : return CoefficientsIndexIterator(l_max_, 0, 0);
294 : }
295 1 : CoefficientsMetadata::CoefficientsIndexIterator cbegin() const {
296 : return begin();
297 : }
298 : /// @}
299 :
300 : /// @{
301 : /// \brief Get a bidirectional iterator to the end of the series of modes.
302 1 : CoefficientsMetadata::CoefficientsIndexIterator end() const {
303 : return CoefficientsIndexIterator(l_max_, l_max_ + 1, l_max_ + 1);
304 : }
305 1 : CoefficientsMetadata::CoefficientsIndexIterator cend() const { return end(); }
306 : /// @}
307 :
308 : private:
309 0 : std::unique_ptr<sharp_alm_info, detail::DestroySharpAlm> alm_info_;
310 0 : size_t l_max_ = 0;
311 : };
312 :
313 : /*!
314 : * \ingroup SwshGroup
315 : * \brief Generation function for obtaining a `CoefficientsMetadata` object
316 : * which is computed by the libsharp calls only once, then lazily cached as a
317 : * singleton via a static member of a function template. This is the preferred
318 : * method for obtaining a `CoefficientsMetadata` when the `l_max` is not very
319 : * large
320 : *
321 : * See the comments in the similar implementation found in `SwshCollocation.hpp`
322 : * for more details on the lazy cache.
323 : */
324 1 : const CoefficientsMetadata& cached_coefficients_metadata(size_t l_max);
325 :
326 : /// @{
327 : /*!
328 : * \ingroup SwshGroup
329 : * \brief Compute the mode coefficient for the convention of
330 : * \cite Goldberg1966uu.
331 : *
332 : * \details See the documentation for `TransformJob` for complete
333 : * details on the libsharp and Goldberg coefficient representations.
334 : * This interface takes a `LibsharpCoefficientInfo`, so that iterating over the
335 : * libsharp data structure and performing computations based on the
336 : * corresponding Goldberg modes is convenient. Two functions are provided, one
337 : * for the positive m modes in the Goldberg representation, and one for the
338 : * negative m modes, as the distinction is not made by the coefficient iterator.
339 : *
340 : * \param coefficient_info An iterator which stores an \f$l\f$ and \f$|m|\f$
341 : * for the mode to be converted to the Goldberg form.
342 : * \param libsharp_modes The libsharp-compatible modal storage to be accessed
343 : * \param radial_offset The index of which angular slice is desired. Modes for
344 : * concentric spherical shells are stored as consecutive blocks of angular
345 : * modes.
346 : */
347 : template <int Spin>
348 1 : std::complex<double> libsharp_mode_to_goldberg_plus_m(
349 : const LibsharpCoefficientInfo& coefficient_info,
350 : const SpinWeighted<ComplexModalVector, Spin>& libsharp_modes,
351 : size_t radial_offset);
352 :
353 : template <int Spin>
354 1 : std::complex<double> libsharp_mode_to_goldberg_minus_m(
355 : const LibsharpCoefficientInfo& coefficient_info,
356 : const SpinWeighted<ComplexModalVector, Spin>& libsharp_modes,
357 : size_t radial_offset);
358 : /// @}
359 :
360 : /*!
361 : * \ingroup SwshGroup
362 : * \brief Compute the mode coefficient for the convention of
363 : * \cite Goldberg1966uu. See the documentation for `TransformJob` for complete
364 : * details on the libsharp and Goldberg coefficient representations.
365 : *
366 : * \details This interface extracts the (`l`, `m`) Goldberg mode from the
367 : * provided libsharp-compatible `libsharp_modes`, at angular slice
368 : * `radial_offset` (Modes for concentric spherical shells are stored as
369 : * consecutive blocks of angular modes). `l_max` must also be specified to
370 : * determine the representation details.
371 : */
372 : template <int Spin>
373 1 : std::complex<double> libsharp_mode_to_goldberg(
374 : size_t l, int m, size_t l_max,
375 : const SpinWeighted<ComplexModalVector, Spin>& libsharp_modes,
376 : size_t radial_offset);
377 :
378 : /*!
379 : * \ingroup SwshGroup
380 : * \brief Set modes of a libsharp-compatible data structure by specifying the
381 : * modes in the \cite Goldberg1966uu representation.
382 : *
383 : * \details This interface takes a `LibsharpCoefficientInfo`, so that iterating
384 : * over the libsharp data structure and performing edits based on the
385 : * corresponding Goldberg modes is convenient.
386 : *
387 : * \param coefficient_info An iterator which stores an \f$l\f$ and \f$|m|\f$
388 : * for the mode to be written by the input Goldberg modes. Note that both the
389 : * \f$+m\f$ and \f$-m\f$ modes must be simultaneously changed, due to the
390 : * representation details. See discussion in `TransformJob` for details.
391 : * \param libsharp_modes The libsharp-compatible modal storage to be altered.
392 : * \param radial_offset The index of which angular slice is desired. Modes for
393 : * concentric spherical shells are stored as consecutive blocks of angular
394 : * modes.
395 : * \param goldberg_plus_m_mode_value The coefficient in the Goldberg
396 : * representation (\f$l\f$, \f$m\f$)
397 : * \param goldberg_minus_m_mode_value The coefficient in the Goldberg
398 : * representation (\f$l\f$, \f$-m\f$)
399 : */
400 : template <int Spin>
401 1 : void goldberg_modes_to_libsharp_modes_single_pair(
402 : const LibsharpCoefficientInfo& coefficient_info,
403 : gsl::not_null<SpinWeighted<ComplexModalVector, Spin>*> libsharp_modes,
404 : size_t radial_offset, std::complex<double> goldberg_plus_m_mode_value,
405 : std::complex<double> goldberg_minus_m_mode_value);
406 :
407 : /*!
408 : * \ingroup SwshGroup
409 : * \brief Set modes of a libsharp-compatible data structure by specifying the
410 : * modes in the \cite Goldberg1966uu representation.
411 : *
412 : * \details This interface sets the (\f$l\f$, \f$\pm m\f$), modes in a
413 : * libsharp-compatible representation from modes specified in the Goldberg
414 : * representation. Note that both the \f$+m\f$ and \f$-m\f$ modes must be
415 : * simultaneously changed, due to the representation details. See discussion in
416 : * `TransformJob` for details.
417 : * \param libsharp_modes The libsharp-compatible modal storage to be altered.
418 : * \param radial_offset The index of which angular slice is desired. Modes for
419 : * concentric spherical shells are stored as consecutive blocks of angular
420 : * modes.
421 : * \param goldberg_plus_m_mode_value The coefficient in the Goldberg
422 : * representation (\f$l\f$, \f$m\f$)
423 : * \param goldberg_minus_m_mode_value The coefficient in the Goldberg
424 : * representation (\f$l\f$, \f$-m\f$)
425 : * \param l the spherical harmonic \f$l\f$ mode requested
426 : * \param m the spherical harmonic \f$m\f$ mode requested
427 : * \param l_max the `l_max` for the provided coefficient set
428 : */
429 : template <int Spin>
430 1 : void goldberg_modes_to_libsharp_modes_single_pair(
431 : size_t l, int m, size_t l_max,
432 : gsl::not_null<SpinWeighted<ComplexModalVector, Spin>*> libsharp_modes,
433 : size_t radial_offset, std::complex<double> goldberg_plus_m_mode_value,
434 : std::complex<double> goldberg_minus_m_mode_value);
435 :
436 : /// @{
437 : /*!
438 : * \ingroup SwshGroup
439 : * \brief Compute the set of Goldberg Spin-weighted spherical harmonic modes (in
440 : * the convention of \cite Goldberg1966uu) from a libsharp-compatible series of
441 : * modes.
442 : *
443 : * \details Modes for concentric spherical shells are stored as consecutive
444 : * blocks of angular modes in both representations.
445 : */
446 : template <int Spin>
447 1 : void libsharp_to_goldberg_modes(
448 : gsl::not_null<SpinWeighted<ComplexModalVector, Spin>*> goldberg_modes,
449 : const SpinWeighted<ComplexModalVector, Spin>& libsharp_modes, size_t l_max);
450 :
451 : template <int Spin>
452 1 : SpinWeighted<ComplexModalVector, Spin> libsharp_to_goldberg_modes(
453 : const SpinWeighted<ComplexModalVector, Spin>& libsharp_modes, size_t l_max);
454 : /// @}
455 :
456 : /// @{
457 : /*!
458 : * \ingroup SwshGroup
459 : * \brief Compute the set of libsharp-compatible spin-weighted spherical
460 : * harmonic modes from a set of Goldberg modes (following the convention of
461 : * \cite Goldberg1966uu)
462 : *
463 : * \details Internally iterates and uses the
464 : * `goldberg_modes_to_libsharp_modes_single_pair()`. In many applications where
465 : * it is not necessary to store the full Goldberg representation, memory
466 : * allocation can be saved by manually performing the iteration and calls to
467 : * `goldberg_modes_to_libsharp_modes_single_pair()`.
468 : */
469 : template <int Spin>
470 1 : void goldberg_to_libsharp_modes(
471 : gsl::not_null<SpinWeighted<ComplexModalVector, Spin>*> libsharp_modes,
472 : const SpinWeighted<ComplexModalVector, Spin>& goldberg_modes, size_t l_max);
473 :
474 : template <int Spin>
475 1 : SpinWeighted<ComplexModalVector, Spin> goldberg_to_libsharp_modes(
476 : const SpinWeighted<ComplexModalVector, Spin>& goldberg_modes, size_t l_max);
477 : /// @}
478 :
479 : /*!
480 : * \ingroup SwshGroup
481 : * \brief Returns the index into a vector of modes consistent with
482 : * \cite Goldberg1966uu.
483 : *
484 : * \details `radial_offset` is used to index into three-dimensional data in
485 : * which concentric spherical shells are stored as consecutive blocks of angular
486 : * modes. Goldberg modes are stored in an m varies fastest, and ascending
487 : * \f$m\f$ and \f$l\f$ values. So, the first several modes are (l, m):
488 : *
489 : * [(0, 0), (1, -1), (1, 0), (1, 1), (2, -2), ...]
490 : */
491 : constexpr SPECTRE_ALWAYS_INLINE size_t
492 1 : goldberg_mode_index(const size_t l_max, const size_t l, const int m,
493 : const size_t radial_offset = 0) {
494 : return static_cast<size_t>(
495 : static_cast<int>(square(l_max + 1) * radial_offset + square(l) + l) + m);
496 : }
497 :
498 : } // namespace Swsh
499 : } // namespace Spectral
|
