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