SwshCoefficients.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/SwshCollocation.hpp"
12 
13 namespace Spectral {
14 namespace Swsh {
15 
16 /// \ingroup SwshGroup
17 /// \brief Convenience function for determining the number of spin-weighted
18 /// spherical harmonics coefficients that are stored for a given `l_max`
19 constexpr SPECTRE_ALWAYS_INLINE size_t
20 number_of_swsh_coefficients(const size_t l_max) noexcept {
21  return (l_max + 1) * (l_max + 2) / 2; // "triangular" representation
22 }
23 
24 /*!
25  * \ingroup SwshGroup
26  * \brief Compute the relative sign change necessary to convert between the
27  * libsharp basis for spin weight `from_spin_weight` to the basis for spin
28  * weight `to_spin_weight`, for the real component coefficients if `real` is
29  * true, otherwise for the imaginary component coefficients. The sign change for
30  * a given coefficient is equivalent to the product of
31  * `sharp_swsh_sign(from_spin, m, real) * sharp_swsh_sign(to_spin, m,
32  * real)`. Due to the form of the signs, it does not end up depending on m (the
33  * m's in the power of \f$-1\f$'s cancel).
34  * For full details of the libsharp sign conventions, see the documentation for
35  * TransformJob.
36  *
37  * \details The sign change is obtained by the
38  * difference between the libsharp convention and the convention which uses:
39  *
40  * \f[
41  * {}_s Y_{\ell m} (\theta, \phi) = (-1)^m \sqrt{\frac{2 l + 1}{4 \pi}}
42  * D^{\ell}{}_{-m s}(\phi, \theta, 0).
43  * \f]
44  *
45  * See \cite Goldberg1966uu
46  */
48  const int from_spin_weight, const int to_spin_weight,
49  const bool real) noexcept {
50  if (real) {
51  return (from_spin_weight == 0 ? -1.0 : 1.0) *
52  (from_spin_weight >= 0 ? -1.0 : 1.0) *
53  ((from_spin_weight < 0 and (from_spin_weight % 2 == 0)) ? -1.0
54  : 1.0) *
55  (to_spin_weight == 0 ? -1.0 : 1.0) *
56  (to_spin_weight >= 0 ? -1.0 : 1.0) *
57  ((to_spin_weight < 0 and (to_spin_weight % 2 == 0)) ? -1.0 : 1.0);
58  }
59  return (from_spin_weight == 0 ? -1.0 : 1.0) *
60  ((from_spin_weight < 0 and (from_spin_weight % 2 == 0)) ? 1.0 : -1.0) *
61  (to_spin_weight == 0 ? -1.0 : 1.0) *
62  ((to_spin_weight < 0 and (to_spin_weight % 2 == 0)) ? 1.0 : -1.0);
63 }
64 
65 /*!
66  * \ingroup SwshGroup
67  * \brief Compute the sign change between the libsharp convention and the set
68  * of spin-weighted spherical harmonics given by the relation to the Wigner
69  * rotation matrices.
70  *
71  * \details The sign change is obtained via the difference between the libsharp
72  * convention and the convention which uses:
73  *
74  * \f[
75  * {}_s Y_{\ell m} (\theta, \phi) = (-1)^m \sqrt{\frac{2 l + 1}{4 \pi}}
76  * D^{\ell}{}_{-m s}(\phi, \theta, 0).
77  * \f]
78  *
79  * See \cite Goldberg1966uu.
80  * The sign change is computed for the real component coefficients if `real` is
81  * true, otherwise for the imaginary component coefficients. For full details on
82  * the sign convention used in libsharp, see the documentation for TransformJob.
83  * This function outputs the \f$\mathrm{sign}(m, s, \mathrm{real})\f$ necessary
84  * to produce the conversion between Goldberg moments and libsharp moments:
85  *
86  * \f{align*}{
87  * {}_s Y_{\ell m}^{\mathrm{real}, \mathrm{sharp}} =& \mathrm{sign}(m, s,
88  * \mathrm{real=true}) {}_s Y_{\ell m}^{\mathrm{Goldberg}}\\
89  * {}_s Y_{\ell m}^{\mathrm{imag}, \mathrm{sharp}} =&
90  * \mathrm{sign}(m, s, \mathrm{real=false}) {}_s Y_{\ell
91  * m}^{\mathrm{Goldberg}}.
92  * \f}
93  *
94  * Note that in this equation, the "real" and "imag" superscripts refer to the
95  * set of basis functions used for the decomposition of the real and imaginary
96  * part of the spin-weighted collocation points, not real or imaginary parts of
97  * the basis functions themselves.
98  */
100  const int spin_weight, const int m, const bool real) noexcept {
101  if (real) {
102  if (m >= 0) {
103  return (spin_weight > 0 ? -1.0 : 1.0) *
104  ((spin_weight < 0 and (spin_weight % 2 == 0)) ? -1.0 : 1.0);
105  }
106  return (spin_weight == 0 ? -1.0 : 1.0) *
107  ((spin_weight >= 0 and (m % 2 == 0)) ? -1.0 : 1.0) *
108  (spin_weight < 0 and (m + spin_weight) % 2 == 0 ? -1.0 : 1.0);
109  }
110  if (m >= 0) {
111  return (spin_weight == 0 ? -1.0 : 1.0) *
112  ((spin_weight < 0 and (spin_weight % 2 == 0)) ? 1.0 : -1.0);
113  }
114  return (spin_weight == 0 ? -1.0 : 1.0) *
115  ((spin_weight >= 0 and (m % 2 == 0)) ? -1.0 : 1.0) *
116  ((spin_weight < 0 and (m + spin_weight) % 2 != 0) ? -1.0 : 1.0);
117 }
118 
119 namespace detail {
120 
121 constexpr size_t coefficients_maximum_l_max = collocation_maximum_l_max;
122 
123 struct DestroySharpAlm {
124  void operator()(sharp_alm_info* to_delete) noexcept {
125  sharp_destroy_alm_info(to_delete);
126  }
127 };
128 // The Coefficients class acts largely as a memory-safe container for a
129 // `sharp_alm_info*`, required for use of libsharp transform utilities.
130 // The libsharp utilities are currently constructed to only provide user
131 // functions with collocation data for spin-weighted functions and
132 // derivatives. If and when the libsharp utilities are expanded to provide
133 // spin-weighted coefficients as output, this class should be expanded to
134 // provide information about the value and storage ordering of those
135 // coefficients to user code. This should be implemented as an iterator, as is
136 // done in SwshCollocation.hpp.
137 //
138 // Note: The libsharp representation of coefficients is altered from the
139 // standard mathematical definitions in a nontrivial way. There are a number of
140 // important features to the data storage of the coefficients.
141 // - they are stored as a set of complex values, but each vector of complex
142 // values is the transform of only the real or imaginary part of the
143 // collocation data
144 // - because each vector of complex coefficients is related to the transform of
145 // a set of doubles, only (about) half of the m's are stored (m >= 0), because
146 // the remaining m modes are determinable by conjugation from the positive m
147 // modes, given that they represent the transform of a purely real or purely
148 // imaginary collocation quantity
149 // - they are stored in an l-varies-fastest, triangular representation. To be
150 // concrete, for an l_max=2, the order of coefficient storage is (l, m):
151 // [(0, 0), (1, 0), (2, 0), (1, 1), (2, 1), (2, 2)]
152 // - due to the restriction of representing only the transform of real
153 // quantities, the m=0 modes always have vanishing imaginary component.
154 class Coefficients {
155  public:
156  explicit Coefficients(size_t l_max) noexcept;
157 
158  ~Coefficients() = default;
159  Coefficients() = default;
160  Coefficients(const Coefficients&) = delete;
161  Coefficients(Coefficients&&) = default;
162  Coefficients& operator=(const Coefficients&) = delete;
163  Coefficients& operator=(Coefficients&&) = default;
164  sharp_alm_info* get_sharp_alm_info() const noexcept {
165  return alm_info_.get();
166  }
167 
168  size_t l_max() const noexcept { return l_max_; }
169 
170  private:
172  size_t l_max_ = 0;
173 };
174 
175 // Function for obtaining a `Coefficients`, which is a thin wrapper around
176 // the libsharp `alm_info`, needed to perform transformations and iterate over
177 // coefficients. A lazy static cache is used to avoid repeated computation. See
178 // the similar implementation in `SwshCollocation.hpp` for details about the
179 // caching mechanism.
180 const Coefficients& precomputed_coefficients(size_t l_max) noexcept;
181 } // namespace detail
182 } // namespace Swsh
183 } // namespace Spectral
constexpr double sharp_swsh_sign_change(const int from_spin_weight, const int to_spin_weight, const bool real) noexcept
Compute the relative sign change necessary to convert between the libsharp basis for spin weight from...
Definition: SwshCoefficients.hpp:47
Definition: Determinant.hpp:11
#define SPECTRE_ALWAYS_INLINE
Always inline a function. Only use this if you benchmarked the code.
Definition: ForceInline.hpp:20
Defines macro to always inline a function.
constexpr size_t number_of_swsh_coefficients(const size_t l_max) noexcept
Convenience function for determining the number of spin-weighted spherical harmonics coefficients tha...
Definition: SwshCoefficients.hpp:20
constexpr double sharp_swsh_sign(const int spin_weight, const int m, const bool real) noexcept
Compute the sign change between the libsharp convention and the set of spin-weighted spherical harmon...
Definition: SwshCoefficients.hpp:99
Functionality associated with a particular choice of basis functions and quadrature for spectral oper...
Definition: Chebyshev.cpp:19