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 <cstddef>
8 : #include <vector>
9 :
10 : #include "DataStructures/ComplexModalVector.hpp"
11 : #include "DataStructures/DataVector.hpp"
12 : #include "DataStructures/SpinWeighted.hpp"
13 : #include "Utilities/Gsl.hpp"
14 :
15 : /// \cond
16 : class ComplexDataVector;
17 : namespace PUP {
18 : class er;
19 : } // namespace PUP
20 : /// \endcond
21 :
22 : namespace Spectral {
23 : namespace Swsh {
24 :
25 : /// \ingroup SwshGroup
26 : /// A utility for evaluating a particular spin-weighted spherical harmonic
27 : /// function at arbitrary points.
28 : ///
29 : /// \warning This should NOT be used for interpolation; such an evaluation
30 : /// strategy is hopelessly slow compared to Clenshaw recurrence strategies.
31 : /// Instead use `SwshInterpolator`.
32 1 : class SpinWeightedSphericalHarmonic {
33 : public:
34 : // charm needs an empty constructor
35 0 : SpinWeightedSphericalHarmonic() = default;
36 :
37 0 : SpinWeightedSphericalHarmonic(int spin, size_t l, int m);
38 :
39 : /*!
40 : * Return by pointer the values of the spin-weighted spherical harmonic
41 : * evaluated at `theta` and `phi`.
42 : *
43 : * \details The additional values `sin_theta_over_2` and `cos_theta_over_2`,
44 : * representing \f$\sin(\theta/2)\f$ and \f$\cos(\theta/2)\f$ are taken as
45 : * required input to improve the speed of evaluation when called more than
46 : * once.
47 : *
48 : * The formula we evaluate (with various prefactors precomputed, cached, and
49 : * optimized from the factorials) is \cite Goldberg1966uu
50 : *
51 : * \f{align*}{
52 : * {}_s Y_{l m} = (-1)^m \sqrt{\frac{(l + m)! (l-m)! (2l + 1)}
53 : * {4 \pi (l + s)! (l - s)!}} \sin^{2 l}(\theta / 2)
54 : * \sum_{r = 0}^{l - s} {l - s \choose r} {l + s \choose r + s - m}
55 : * (-1)^{l - r - s} e^{i m \phi} \cot^{2 r + s - m}(\theta / 2).
56 : * \f}
57 : */
58 1 : void evaluate(gsl::not_null<ComplexDataVector*> result,
59 : const DataVector& theta, const DataVector& phi,
60 : const DataVector& sin_theta_over_2,
61 : const DataVector& cos_theta_over_2) const;
62 :
63 : /// Return by value the spin-weighted spherical harmonic evaluated at `theta`
64 : /// and `phi`.
65 : ///
66 : /// \details The additional values `sin_theta_over_2` and `cos_theta_over_2`,
67 : /// representing \f$\sin(\theta/2)\f$ and \f$\cos(\theta/2)\f$ are taken as
68 : /// required input to improve the speed of evaluation when called more than
69 : /// once.
70 1 : ComplexDataVector evaluate(const DataVector& theta, const DataVector& phi,
71 : const DataVector& sin_theta_over_2,
72 : const DataVector& cos_theta_over_2) const;
73 :
74 : /// Return by value the spin-weighted spherical harmonic evaluated at `theta`
75 : /// and `phi`.
76 1 : std::complex<double> evaluate(double theta, double phi) const;
77 :
78 : /// Serialization for Charm++.
79 : // NOLINTNEXTLINE(google-runtime-references)
80 1 : void pup(PUP::er& p);
81 :
82 : private:
83 0 : int spin_ = 0;
84 0 : size_t l_ = 0;
85 0 : int m_ = 0;
86 0 : double overall_prefactor_ = 0.0;
87 0 : std::vector<double> r_prefactors_;
88 : };
89 :
90 : /*!
91 : * \ingroup SwshGroup
92 : * \brief Performs interpolation for spin-weighted spherical harmonics by
93 : * taking advantage of the Clenshaw method of expanding recurrence relations.
94 : *
95 : * \details During construction, we cache several functions of the target
96 : * interpolation points that will be used during the Clenshaw evaluation.
97 : * A new `SwshInterpolator` object must be created for each new set of target
98 : * points, but the member function `SwshInterpolator::interpolate()` may be
99 : * called on several different coefficients or collocation sets, and of
100 : * different spin-weights.
101 : *
102 : * Recurrence constants
103 : * --------------------
104 : * This utility obtains the Clenshaw interpolation constants from a
105 : * `StaticCache`, so that 'universal' quantities can be calculated only once per
106 : * execution and re-used on each interpolation.
107 : *
108 : * We evaluate the recurrence coefficients \f$\alpha_l^{(a,b)}\f$ and
109 : * \f$\beta_l^{(a,b)}\f$, where \f$a = |s + m|\f$, \f$b = |s - m|\f$, and
110 : *
111 : * \f[
112 : * {}_sY_{l, m}(\theta, \phi) = \alpha_l^{(a,b)}(\theta) {}_s Y_{l-1, m}
113 : * + \beta_l^{(a, b)} {}_s Y_{l - 2, m}(\theta, \phi)
114 : * \f]
115 : *
116 : * The core Clenshaw recurrence is in the\f$l\f$ modes of the
117 : * spin-weighted spherical harmonic. For a set of modes \f$a_l^{(a,b)}\f$ at
118 : * fixed \f$m\f$, the function value is evaluated by the recurrence:
119 : *
120 : * \f{align*}{
121 : * y^{(a, b)}_{l_\text{max} + 2} &= y^{(a, b)}_{l_\text{max} + 1} = 0 \\
122 : * y^{(a, b)}_{l}(\theta) &= \alpha_{l + 1}^{(a, b)} y^{(a, b)}_{l +
123 : * 1}(\theta)
124 : * + \beta_{l + 2}^{(a,b)} y^{(a, b)}_{l + 2}(\theta) + a_l^{(a, b)} \\
125 : * f_m(\theta, \phi) &= \beta_{l_{\text{min}} + 2}
126 : * {}_s Y_{l_{\text{min}}, m}(\theta, \phi) y^{(a, b)}_2(\theta) +
127 : * {}_s Y_{l_{\text{min}} + 1, m}(\theta, \phi) y^{(a, b)}_1(\theta) +
128 : * a_0^{(a, b)} {}_s Y_{l_{\text{min}}, m}(\theta, \phi),
129 : * \f}
130 : *
131 : * where \f$l_{\text{min}} = \max(|m|, |s|)\f$ is the lowest nonvanishing
132 : * \f$l\f$ mode for a given \f$m\f$ and \f$s\f$.
133 : *
134 : * The coefficients to cache are inferred from a mechanical but lengthy
135 : * calculation involving the recurrence relation for the Jacobi polynomials.
136 : * The result is:
137 : *
138 : * \f{align*}{
139 : * \alpha_l^{(a, b)} &= \frac{\sqrt{(2l + 1)(2l - 1)}}
140 : * {2\sqrt{(l + k)(l + k + a + b)(l + k + a)(l + k + b)}}
141 : * \left[2 l \cos(\theta) +
142 : * \frac{a^2 - b^2}{2 l - 2}\right] \\
143 : * \beta_l^{(a, b)} & = - \sqrt{\frac{(2 l + 1)(l + k + a - 1)(l + k + b - 1)
144 : * (l + k - 1)(l + k + a + b - 1)}
145 : * {(2l - 3)(l + k)(l + k + a + b)(l + k + a)(l + k + b)}}
146 : * \frac{2 l}{2 l - 2},
147 : * \f}
148 : * where \f$k = - (a + b)/2\f$ (which is always integral due to the properties
149 : * of \f$a\f$ and \f$b\f$). Note that because the values of \f$\alpha\f$ and
150 : * \f$\beta\f$ in the recurrence relation are not needed for any value below
151 : * \f$l_{\text{min}} + 2\f$, so none of the values in the square-roots or
152 : * denominators take pathological values for any of the coefficients we require.
153 : *
154 : * The \f$\beta\f$ constants are filled in member variable `beta_constant`. The
155 : * \f$\alpha\f$ values are separately stored as the prefactor for the
156 : * \f$\cos(\theta)\f$ term and the constant term in `alpha_prefactor` and
157 : * `alpha_constant` member variables.
158 : *
159 : * In addition, it is efficient to cache recurrence coefficients necessary for
160 : * generating the first couple of spin-weighted spherical harmonic functions for
161 : * each \f$m\f$ used in the Clenshaw sum.
162 : *
163 : * The member variable `harmonic_at_l_min_prefactors` holds the prefactors for
164 : * directly evaluating the harmonics at \f$s >= m\f$,
165 : *
166 : * \f{align*}{
167 : * {}_s Y_{|s| m} = {}_s C_{m} e^{i m \phi} \sin^a(\theta/2) \cos^b(\theta/2),
168 : * \f}
169 : *
170 : * where \f${}_s C_m\f$ are the cached prefactors that take the values
171 : *
172 : * \f{align*}{
173 : * {}_s C_m = (-1)^{m + \lambda(m)} \sqrt{\frac{(2 |s| + 1) (|s| + k)!}
174 : * {4 \pi (|s| + k + a)! (|s| + k + b)!}}
175 : * \f}
176 : *
177 : * and
178 : *
179 : * \f{align*}{
180 : * \lambda(m) =
181 : * \left\{
182 : * \begin{array}{ll}
183 : * 0 &\text{ for } s \ge -m\\
184 : * s + m &\text{ for } s < -m
185 : * \end{array} \right..
186 : * \f}
187 : *
188 : * The member variable `harmonic_m_recurrence_prefactors` holds the prefactors
189 : * necessary to evaluate the lowest harmonics for each \f$m\f$ from the next
190 : * lowest-in-magnitude \f$m\f$, allowing most leading harmonics to recursively
191 : * derived.
192 : *
193 : * \f{align*}{
194 : * {}_s Y_{|m|, m} = {}_s D_m \sin(\theta / 2) \cos(\theta / 2)
195 : * \left\{
196 : * \begin{array}{ll}
197 : * e^{i \phi} {}_s Y_{|m| - 1, m - 1} &\text{ for } m > 0\\
198 : * e^{-i \phi} {}_s Y_{|m| - 1, m + 1} &\text{ for } m < 0
199 : * \end{array}\right.,
200 : * \f}
201 : *
202 : * where \f${}_s D_m\f$ are the cached prefactors that take the values
203 : *
204 : * \f{align*}{
205 : * {}_s D_m = -(-1)^{\Delta \lambda(m)} \sqrt{
206 : * \frac{(2 |m| + 1)(k + |m| + a + b - 1)(k + |m| + a + b)}
207 : * {(2 |m| - 1)(k + |m| + a)(k + |m| + b)}},
208 : * \f}
209 : * and \f$\Delta \lambda(m)\f$ is the difference in \f$\lambda(m)\f$ between
210 : * the harmonic on the left-hand side and right-hand side of the above
211 : * recurrence relation, that is \f$\lambda(m) - \lambda(m - 1)\f$ for positive
212 : * \f$m\f$ and \f$\lambda(m) - \lambda(m + 1)\f$ for negative \f$m\f$.
213 : *
214 : * Finally, the member variable
215 : * `harmonic_at_l_min_plus_one_recurrence_prefactors` holds the prefactors
216 : * necessary to evaluate parts of the recurrence relations from the lowest
217 : * \f$l\f$ for a given \f$m\f$ to the next-to-lowest \f$l\f$.
218 : *
219 : * \f{align*}{
220 : * {}_s Y_{l_{\min} + 1, m} = {}_s E_m
221 : * \left[(a + 1) + (a + b + 2) \frac{(\cos(\theta) - 1)}{2}\right]
222 : * {}_s Y_{l_{\min}, m}.
223 : * \f}
224 : *
225 : * where \f${}_s E_m\f$ are the cached prefactors that take the values
226 : *
227 : * \f{align*}{
228 : * {}_s E_{m} = \sqrt{\frac{2 l_{\min} + 3}{ 2 l_{\min} + 1}}
229 : * \sqrt{\frac{(l_{\min} + k + 1)(l_{\min} + k + a + b + 1)}
230 : * {(l_{\min} + k + a + 1)(l_{\min} + k + b + 1)}}
231 : * \f}
232 : */
233 1 : class SwshInterpolator {
234 : public:
235 : // charm needs the empty constructor
236 0 : SwshInterpolator() = default;
237 :
238 0 : SwshInterpolator(const SwshInterpolator&) = default;
239 0 : SwshInterpolator(SwshInterpolator&&) = default;
240 0 : SwshInterpolator& operator=(const SwshInterpolator&) = default;
241 0 : SwshInterpolator& operator=(SwshInterpolator&&) = default;
242 0 : ~SwshInterpolator() = default;
243 :
244 0 : SwshInterpolator(const DataVector& theta, const DataVector& phi,
245 : size_t l_max);
246 :
247 : /*!
248 : * \brief Perform the Clenshaw recurrence sum, returning by pointer
249 : * `interpolated` of interpolating the `goldberg_modes` at the collocation
250 : * points passed to the constructor.
251 : *
252 : * \details The core Clenshaw recurrence is in the\f$l\f$ modes of the
253 : * spin-weighted spherical harmonic. For a set of modes \f$a_l^{(a,b)}\f$ at
254 : * fixed \f$m\f$, the function value is evaluated by the recurrence:
255 : *
256 : * \f{align*}{
257 : * y^{(a, b)}_{l_\text{max} + 2} &= y^{(a, b)}_{l_\text{max} + 1} = 0 \\
258 : * y^{(a, b)}_{l}(\theta) &= \alpha_{l + 1}^{(a, b)} y^{(a, b)}_{l +
259 : * 1}(\theta)
260 : * + \beta_{l + 2}^{(a,b)} y^{(a, b)}_{l + 2}(\theta) + a_l^{(a, b)} \\
261 : * f_m(\theta, \phi) &= \beta_{l_{\text{min}} + 2}
262 : * {}_s Y_{l_{\text{min}}, m}(\theta, \phi) y^{(a, b)}_2(\theta) +
263 : * {}_s Y_{l_{\text{min}} + 1, m}(\theta, \phi) y^{(a, b)}_1(\theta) +
264 : * a_0^{(a, b)} {}_s Y_{l_{\text{min}}, m}(\theta, \phi)
265 : * \f}
266 : *
267 : * The recurrence in \f$l\f$ accomplishes much of the work, but for full
268 : * efficiency, we also recursively evaluate the lowest handful of \f$l\f$s for
269 : * each \f$m\f$. The details of those additional recurrence tricks can be
270 : * found in the documentation for `ClenshawRecurrenceConstants`.
271 : */
272 : template <int Spin>
273 1 : void interpolate(
274 : gsl::not_null<SpinWeighted<ComplexDataVector, Spin>*> interpolated,
275 : const SpinWeighted<ComplexModalVector, Spin>& goldberg_modes) const;
276 :
277 : /*!
278 : * \brief Perform the Clenshaw recurrence sum, returning by pointer
279 : * `interpolated` of interpolating function represented by
280 : * `libsharp_collocation` at the target points passed to the constructor.
281 : *
282 : * \details The core Clenshaw recurrence is in the\f$l\f$ modes of the
283 : * spin-weighted spherical harmonic. For a set of modes \f$a_l^{(a,b)}\f$ at
284 : * fixed \f$m\f$, the function value is evaluated by the recurrence:
285 : *
286 : * \f{align*}{
287 : * y^{(a, b)}_{l_\text{max} + 2} &= y^{(a, b)}_{l_\text{max} + 1} = 0 \\
288 : * y^{(a, b)}_{l}(\theta) &= \alpha_{l + 1}^{(a, b)} y^{(a, b)}_{l +
289 : * 1}(\theta)
290 : * + \beta_{l + 2}^{(a,b)} y^{(a, b)}_{l + 2}(\theta) + a_l^{(a, b)} \\
291 : * f_m(\theta, \phi) &= \beta_{l_{\text{min}} + 2}
292 : * {}_s Y_{l_{\text{min}}, m}(\theta, \phi) y^{(a, b)}_2(\theta) +
293 : * {}_s Y_{l_{\text{min}} + 1, m}(\theta, \phi) y^{(a, b)}_1(\theta) +
294 : * a_0^{(a, b)} {}_s Y_{l_{\text{min}}, m}(\theta, \phi)
295 : * \f}
296 : *
297 : * The recurrence in \f$l\f$ accomplishes much of the work, but for full
298 : * efficiency, we also recursively evaluate the lowest handful of \f$l\f$s for
299 : * each \f$m\f$. The details of those additional recurrence tricks can be
300 : * found in the documentation for `ClenshawRecurrenceConstants`.
301 : */
302 : template <int Spin>
303 1 : void interpolate(
304 : gsl::not_null<SpinWeighted<ComplexDataVector, Spin>*> interpolated,
305 : const SpinWeighted<ComplexDataVector, Spin>& libsharp_collocation) const;
306 :
307 : /// \brief Evaluate the SWSH function at the lowest \f$l\f$ value for a given
308 : /// \f$m\f$ at the target interpolation points.
309 : ///
310 : /// \details Included in the public interface for thorough testing, most use
311 : /// cases should just use the `SwshInterpolator::interpolate()` member
312 : /// function.
313 : template <int Spin>
314 1 : void direct_evaluation_swsh_at_l_min(
315 : gsl::not_null<SpinWeighted<ComplexDataVector, Spin>*> harmonic,
316 : int m) const;
317 :
318 : /// \brief Evaluate the SWSH function at the next-to-lowest \f$l\f$ value for
319 : /// a given \f$m\f$ at the target interpolation points, given input harmonic
320 : /// values for the lowest \f$l\f$ value.
321 : ///
322 : /// \details Included in the public interface for thorough testing, most use
323 : /// cases should just use the `SwshInterpolator::interpolate()` member
324 : /// function.
325 : template <int Spin>
326 1 : void evaluate_swsh_at_l_min_plus_one(
327 : gsl::not_null<SpinWeighted<ComplexDataVector, Spin>*> harmonic,
328 : const SpinWeighted<ComplexDataVector, Spin>& harmonic_at_l_min,
329 : int m) const;
330 :
331 : /// \brief Evaluate the SWSH function at the lowest \f$l\f$ value for a given
332 : /// \f$m\f$ at the target interpolation points, given harmonic data at the
333 : /// next lower \f$m\f$ (by magnitude), passed in by the same pointer used for
334 : /// the return.
335 : ///
336 : /// \details Included in the public interface for thorough testing, most use
337 : /// cases should just use the `SwshInterpolator::interpolate()` member
338 : /// function.
339 : template <int Spin>
340 1 : void evaluate_swsh_m_recurrence_at_l_min(
341 : gsl::not_null<SpinWeighted<ComplexDataVector, Spin>*> harmonic,
342 : int m) const;
343 :
344 : /// \brief Perform the core Clenshaw interpolation at fixed \f$m\f$,
345 : /// accumulating the result in `interpolation`.
346 : ///
347 : /// \details Included in the public interface for thorough testing, most use
348 : /// cases should just use the `interpolate` member function.
349 : template <int Spin>
350 1 : void clenshaw_sum(
351 : gsl::not_null<SpinWeighted<ComplexDataVector, Spin>*> interpolation,
352 : const SpinWeighted<ComplexDataVector, Spin>& l_min_harmonic,
353 : const SpinWeighted<ComplexDataVector, Spin>& l_min_plus_one_harmonic,
354 : const SpinWeighted<ComplexModalVector, Spin>& goldberg_modes,
355 : int m) const;
356 :
357 : /// Serialization for Charm++.
358 : // NOLINTNEXTLINE(google-runtime-references)
359 1 : void pup(PUP::er& p);
360 :
361 : private:
362 0 : size_t l_max_ = 0;
363 0 : DataVector cos_theta_;
364 0 : DataVector sin_theta_;
365 0 : DataVector cos_theta_over_two_;
366 0 : DataVector sin_theta_over_two_;
367 0 : std::vector<DataVector> sin_m_phi_;
368 0 : std::vector<DataVector> cos_m_phi_;
369 : // NOLINTNEXTLINE(spectre-mutable)
370 0 : mutable ComplexModalVector raw_libsharp_coefficient_buffer_;
371 : // NOLINTNEXTLINE(spectre-mutable)
372 0 : mutable ComplexModalVector raw_goldberg_coefficient_buffer_;
373 : };
374 : } // namespace Swsh
375 : } // namespace Spectral
|
