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