SwshTransform.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 <sharp_cxx.h>
9 #include <type_traits>
10 #include <utility>
11 #include <vector>
12 
14 #include "DataStructures/SpinWeighted.hpp"
15 #include "DataStructures/Tags.hpp"
16 #include "DataStructures/Variables.hpp" // IWYU pragma: keep
17 #include "NumericalAlgorithms/Spectral/ComplexDataView.hpp"
18 #include "NumericalAlgorithms/Spectral/SwshCoefficients.hpp" // IWYU pragma: keep
19 #include "NumericalAlgorithms/Spectral/SwshCollocation.hpp"
20 #include "NumericalAlgorithms/Spectral/SwshTags.hpp"
21 #include "Utilities/Gsl.hpp"
22 #include "Utilities/TMPL.hpp"
23 
24 // IWYU pragma: no_forward_declare Variables
25 // IWYU pragma: no_forward_declare Coefficients
26 // IWYU pragma: no_forward_declare SpinWeighted
27 
28 class ComplexDataVector;
29 class ComplexModalVector;
30 
31 namespace Spectral {
32 namespace Swsh {
33 
34 namespace detail {
35 // libsharp has an internal maximum number of transforms that is not in the
36 // public interface, so we must hard-code its value here
37 static const size_t max_libsharp_transforms = 100;
38 
39 // appends to a vector of pointers `collocation_data` the set of pointers
40 // associated with angular views of the provided collocation data `vector`. The
41 // associated views are appended into `collocation_views`. If `positive_spin` is
42 // false, this function will conjugate the views, and therefore potentially
43 // conjugate (depending on `Representation`) the input data `vectors`. This
44 // behavior is chosen to avoid frequent copies of the input data `vector` for
45 // optimization. The conjugation must be undone after the call to libsharp
46 // transforms using `conjugate_views` to restore the input data to its original
47 // state.
48 template <ComplexRepresentation Representation>
49 void append_libsharp_collocation_pointers(
50  gsl::not_null<std::vector<double*>*> collocation_data,
51  gsl::not_null<std::vector<ComplexDataView<Representation>>*>
52  collocation_views,
53  gsl::not_null<ComplexDataVector*> vector, size_t l_max,
54  bool positive_spin) noexcept;
55 
56 // Perform a conjugation on a vector of `ComplexDataView`s if `Spin` < 0. This
57 // is used for undoing the conjugation described in the code comment for
58 // `append_libsharp_collocation_pointers`
59 template <int Spin, ComplexRepresentation Representation>
60 SPECTRE_ALWAYS_INLINE void conjugate_views(
61  gsl::not_null<std::vector<ComplexDataView<Representation>>*>
62  collocation_views) noexcept {
63  if (Spin < 0) {
64  for (auto& view : *collocation_views) {
65  view.conjugate();
66  }
67  }
68 }
69 
70 // appends to a vector of pointers the set of pointers associated with
71 // angular views of the provided coefficient vector. When working with the
72 // libsharp coefficient representation, note the intricacies mentioned in
73 // the documentation for `TransformJob`.
74 void append_libsharp_coefficient_pointers(
75  gsl::not_null<std::vector<std::complex<double>*>*> coefficient_data,
76  gsl::not_null<ComplexModalVector*> vector, size_t l_max) noexcept;
77 
78 // perform the actual libsharp execution calls on an input and output set of
79 // pointers. This function handles the complication of a limited maximum number
80 // of simultaneous transforms, performing multiple execution calls on pointer
81 // blocks if necessary.
82 template <ComplexRepresentation Representation>
83 void execute_libsharp_transform_set(
84  const sharp_jobtype& jobtype, int spin,
85  gsl::not_null<std::vector<std::complex<double>*>*> coefficient_data,
86  gsl::not_null<std::vector<double*>*> collocation_data,
87  gsl::not_null<const CollocationMetadata<Representation>*>
88  collocation_metadata,
89  const sharp_alm_info* alm_info, size_t num_transforms) noexcept;
90 
91 // template 'implementation' for the `swsh_transform` function below which
92 // performs an arbitrary number of transforms, and places them in the same
93 // number of destination modal containers passed by pointer.
94 // template notes: the parameter pack represents a collection of ordered modal
95 // data passed by pointer, followed by ordered nodal data passed by const
96 // reference. The first modal value is separated to infer the `Spin` template
97 // parameter.
98 template <ComplexRepresentation Representation, int Spin,
99  typename... ModalThenNodalTypes, size_t... Is>
100 void swsh_transform_impl(
101  size_t l_max, size_t number_of_radial_points,
104  const ModalThenNodalTypes&... coefficients_then_collocations) noexcept;
105 
106 // template 'implementation' for the `inverse_swsh_transform` function below
107 // which performs an arbitrary number of transforms, and places them in the same
108 // number of destination nodal containers passed by pointer.
109 // template notes: the parameter pack represents a collection of ordered nodal
110 // data passed by pointer, followed by ordered modal data passed by const
111 // reference. The first nodal value is separated to infer the `Spin` template
112 // parameter.
113 template <ComplexRepresentation Representation, int Spin,
114  typename... NodalThenModalTypes, size_t... Is>
115 void inverse_swsh_transform_impl(
116  size_t l_max, size_t number_of_radial_points,
119  const NodalThenModalTypes&... collocations_then_coefficients) noexcept;
120 
121 // forward declaration for friend specification of helper from
122 // `SwshDerivatives.hpp`
123 template <typename Transform, typename FullTagList>
124 struct dispatch_to_transform;
125 } // namespace detail
126 
127 /*!
128  * \ingroup SwshGroup
129  * \brief Perform a forward libsharp spin-weighted spherical harmonic transform
130  * on any number of supplied `SpinWeighted<ComplexDataVector, Spin>`.
131  *
132  * \details This function places the result in one or more
133  * `SpinWeighted<ComplexModalVector, Spin>` returned via provided
134  * pointer. This is a simpler interface to the same functionality as the
135  * \ref DataBoxGroup mutation compatible `SwshTransform`.
136  *
137  * The following parameter ordering for the multiple input interface is enforced
138  * by the interior function calls, but is not obvious from the explicit
139  * parameter packs in this function signature:
140  *
141  * - `size_t l_max` : angular resolution for the transform
142  * - `size_t number_of_radial_points` : radial resolution (number of consecutive
143  * transforms in each of the vectors)
144  * - any number of `gsl::not_null<SpinWeighted<ComplexModalVector,
145  * Spin>*>`, all sharing the same `Spin` : The return-by-pointer containers for
146  * the transformed modal quantities
147  * - the same number of `const SpinWeighted<ComplexDataVector, Spin>&`, with the
148  * same `Spin` as the previous function argument set : The input containers of
149  * nodal spin-weighted spherical harmonic data.
150  *
151  * template parameters:
152  * - `Representation`: Either `ComplexRepresentation::Interleaved` or
153  * `ComplexRepresentation::RealsThenImags`, indicating the representation for
154  * intermediate steps of the transformation. The two representations will give
155  * identical results but may help or hurt performance depending on the task.
156  * If unspecified, defaults to `ComplexRepresentation::Interleaved`.
157  *
158  * The result is a set of libsharp-compatible coefficients.
159  * \see SwshTransform for more details on the mathematics of the libsharp
160  * data structures.
161  *
162  * \warning The collocation data is taken by const reference, but can be
163  * temporarily altered in-place during intermediate parts of the computation.
164  * The input data is guaranteed to return to its original state by the end of
165  * the function. In a setting in which multiple threads access the same data
166  * passed as input to this function, a lock must be used to prevent access
167  * during the execution of the transform.
168  */
169 template <
170  ComplexRepresentation Representation = ComplexRepresentation::Interleaved,
171  int Spin, typename... ModalThenNodalTypes>
173  const size_t l_max, const size_t number_of_radial_points,
175  first_coefficient,
176  const ModalThenNodalTypes&... coefficients_then_collocations) noexcept {
177  static_assert(
178  sizeof...(ModalThenNodalTypes) % 2 == 1,
179  "The number of modal arguments passed to swsh_transform must equal the "
180  "number of nodal arguments, so the total number of arguments must be "
181  "even.");
182  detail::swsh_transform_impl<Representation>(
183  l_max, number_of_radial_points,
184  std::make_index_sequence<(sizeof...(coefficients_then_collocations) + 1) /
185  2>{},
186  first_coefficient, coefficients_then_collocations...);
187 }
188 
189 /*!
190  * \ingroup SwshGroup
191  * \brief Perform a forward libsharp spin-weighted spherical harmonic transform
192  * on a single supplied `SpinWeighted<ComplexDataVector, Spin>`.
193  *
194  * \details This function returns a `SpinWeighted<ComplexModalVector, Spin>` by
195  * value (causing an allocation). This is a simpler interface to the same
196  * functionality as the \ref DataBoxGroup mutation compatible `SwshTransform`.
197  * If you have two or more vectors to transform at once, consider the
198  * pass-by-pointer version of `Spectral::Swsh::swsh_transform` or the \ref
199  * DataBoxGroup interface `InverseSwshTransform`, as they are more efficient for
200  * performing several transforms at once.
201  *
202  * The result is a set of libsharp-compatible coefficients.
203  * \see SwshTransform for more details on the mathematics of the libsharp
204  * data structures.
205  *
206  * \warning The collocation data is taken by const reference, but can be
207  * temporarily altered in-place during intermediate parts of the computation.
208  * The input data is guaranteed to return to its original state by the end of
209  * the function. In a setting in which multiple threads access the same data
210  * passed as input to this function, a lock must be used to prevent access
211  * during the execution of the transform.
212  */
213 template <
214  ComplexRepresentation Representation = ComplexRepresentation::Interleaved,
215  int Spin>
217  size_t l_max, size_t number_of_radial_points,
218  const SpinWeighted<ComplexDataVector, Spin>& collocation) noexcept;
219 
220 /*!
221  * \ingroup SwshGroup
222  * \brief Perform an inverse libsharp spin-weighted spherical harmonic transform
223  * on any number of supplied `SpinWeighted<ComplexModalVector, Spin>`.
224  *
225  * \details This function places the result in one or more
226  * `SpinWeighted<ComplexDataVector, Spin>` returned via provided
227  * pointer. This is a simpler interface to the same functionality as the
228  * \ref DataBoxGroup mutation compatible `InverseSwshTransform`.
229  *
230  * The following parameter ordering for the multiple input interface is enforced
231  * by the interior function calls, but is not obvious from the explicit
232  * parameter packs in this function signature:
233  *
234  * - `size_t l_max` : angular resolution for the transform
235  * - `size_t number_of_radial_points` : radial resolution (number of consecutive
236  * transforms in each of the vectors)
237  * - any number of `gsl::not_null<SpinWeighted<ComplexDataVector,
238  * Spin>*>`, all sharing the same `Spin` : The return-by-pointer containers for
239  * the transformed nodal quantities
240  * - the same number of `const SpinWeighted<ComplexModalVector, Spin>&`, with
241  * the same `Spin` as the previous function argument set : The input containers
242  * of modal spin-weighted spherical harmonic data.
243  *
244  * template parameters:
245  * - `Representation`: Either `ComplexRepresentation::Interleaved` or
246  * `ComplexRepresentation::RealsThenImags`, indicating the representation for
247  * intermediate steps of the transformation. The two representations will give
248  * identical results but may help or hurt performance depending on the task.
249  * If unspecified, defaults to `ComplexRepresentation::Interleaved`.
250  *
251  * The input is a set of libsharp-compatible coefficients.
252  * \see SwshTransform for more details on the mathematics of the libsharp
253  * data structures.
254  */
255 template <
256  ComplexRepresentation Representation = ComplexRepresentation::Interleaved,
257  int Spin, typename... NodalThenModalTypes>
259  const size_t l_max, const size_t number_of_radial_points,
261  first_collocation,
262  const NodalThenModalTypes&... collocations_then_coefficients) noexcept {
263  static_assert(
264  sizeof...(NodalThenModalTypes) % 2 == 1,
265  "The number of nodal arguments passed to inverse_swsh_transform must "
266  "equal the number of modal arguments, so the total number of arguments "
267  "must be even.");
268  detail::inverse_swsh_transform_impl<Representation>(
269  l_max, number_of_radial_points,
270  std::make_index_sequence<(sizeof...(collocations_then_coefficients) + 1) /
271  2>{},
272  first_collocation, collocations_then_coefficients...);
273 }
274 
275 /*!
276  * \ingroup SwshGroup
277  * \brief Perform an inverse libsharp spin-weighted spherical harmonic transform
278  * on a single supplied `SpinWeighted<ComplexModalVector, Spin>`.
279  *
280  * \details This function returns a `SpinWeighted<ComplexDataVector, Spin>` by
281  * value (causing an allocation). This is a simpler interface to the same
282  * functionality as the \ref DataBoxGroup mutation compatible
283  * `InverseSwshTransform`. If you have two or more vectors to transform at once,
284  * consider the pass-by-pointer version of
285  * `Spectral::Swsh::inverse_swsh_transform` or the \ref DataBoxGroup interface
286  * `SwshTransform`, as they are more efficient for performing several transforms
287  * at once.
288  *
289  * The input is a set of libsharp-compatible coefficients.
290  * \see SwshTransform for more details on the mathematics of the libsharp
291  * data structures.
292  */
293 template <
294  ComplexRepresentation Representation = ComplexRepresentation::Interleaved,
295  int Spin>
297  size_t l_max, size_t number_of_radial_points,
299  libsharp_coefficients) noexcept;
300 
301 /*!
302  * \ingroup SwshGroup
303  * \brief A \ref DataBoxGroup mutate-compatible computational struct for
304  * performing several spin-weighted spherical harmonic transforms. Internally
305  * dispatches to libsharp.
306  *
307  * \details
308  * Template Parameters:
309  * - `Representation`: The element of the `ComplexRepresentation` enum which
310  * parameterizes how the data passed to libsharp will be represented. Two
311  * options are available:
312  * - `ComplexRepresentation:Interleaved`: indicates that the real and imaginary
313  * parts of the collocation values will be passed to libsharp as pointers
314  * into existing complex data structures, minimizing copies, but maintaining
315  * a stride of 2 for 'adjacent' real or imaginary values.
316  * - `ComplexRepresentation::RealsThenImags`: indicates that the real and
317  * imaginary parts of the collocation values will be passed to libsharp as
318  * separate contiguous blocks. At current, this introduces both allocations
319  * and copies. **optimization note** In the future most of the allocations
320  * can be aggregated to calling code, which would pass in buffers by
321  * `not_null` pointers.
322  *
323  * For performance-sensitive code, both options should be tested, as each
324  * strategy has trade-offs.
325  * - `TagList`: A `tmpl::list` of Tags to be forward transformed. The tags must
326  * represent the nodal data.
327  *
328  * \note The signs obtained from libsharp transformations must be handled
329  * carefully. (In particular, it does not use the sign convention you will find
330  * in [Wikipedia]
331  * (https://en.wikipedia.org/wiki/Spin-weighted_spherical_harmonics)).
332  * - libsharp deals only with the transformation of real values, so
333  * transformation of complex values must be done for real and imaginary parts
334  * separately.
335  * - due to only transforming real components, it stores only the positive
336  * \f$m\f$ modes (as the rest would be redundant). Therefore, the negative
337  * \f$m\f$ modes must be inferred from conjugation and further sign changes.
338  * - libsharp only has the capability of transforming positive spin weighted
339  * quantities. Therefore, additional steps are taken (involving further
340  * conjugation of the provided data) in order to use those utilities on
341  * negative spin-weighted quantities. The resulting modes have yet more sign
342  * changes that must be taken into account.
343  *
344  * The decomposition resulting from the libsharp transform for spin \f$s\f$ and
345  * complex spin-weighted \f${}_s f\f$ can be represented mathematically as:
346  *
347  * \f{align*}{
348  * {}_s f(\theta, \phi) = \sum_{\ell = 0}^{\ell_\mathrm{max}} \Bigg\{&
349  * \left(\sum_{m = 0}^{\ell} a^\mathrm{sharp, real}_{l m} {}_s Y_{\ell
350  * m}^\mathrm{sharp, real}\right) + \left(\sum_{m=1}^{\ell}
351  * \left(a^\mathrm{sharp, real}_{l m}{}\right)^*
352  * {}_s Y_{\ell -m}^\mathrm{sharp, real}\right) \notag\\
353  * &+ i \left[\left(\sum_{m = 0}^{\ell} a^\mathrm{sharp, imag}_{l m} {}_s
354  * Y_{\ell m}^\mathrm{sharp, imag}\right) + \left(\sum_{m=1}^{\ell}
355  * \left(a^\mathrm{sharp, imag}_{l m}{}\right)^* {}_s Y_{\ell -m}^\mathrm{sharp,
356  * imag} \right)\right] \Bigg\},
357  * \f}
358  *
359  * where
360  *
361  * \f{align*}{
362  * {}_s Y_{\ell m}^\mathrm{sharp, real} &=
363  * \begin{cases}
364  * (-1)^{s + 1} {}_s Y_{\ell m}, & \mathrm{for}\; s < 0, m \ge 0 \\
365  * {}_s Y_{\ell m}, & \mathrm{for}\; s = 0, m \ge 0 \\
366  * - {}_s Y_{\ell m}, & \mathrm{for}\; s > 0, m \ge 0 \\
367  * (-1)^{s + m + 1} {}_s Y_{\ell m}, & \mathrm{for}\; s < 0, m < 0 \\
368  * (-1)^{m} {}_s Y_{\ell m}, & \mathrm{for}\; s = 0, m < 0 \\
369  * (-1)^{m + 1} {}_s Y_{\ell m}, & \mathrm{for}\; s > 0, m < 0
370  * \end{cases} \\
371  * &\equiv {}_s S_{\ell m}^{\mathrm{real}} {}_s Y_{\ell m}\\
372  * {}_s Y_{\ell m}^\mathrm{sharp, imag} &=
373  * \begin{cases}
374  * (-1)^{s + 1} {}_s Y_{\ell m}, & \mathrm{for}\; s < 0, m \ge 0 \\
375  * -{}_s Y_{\ell m}, & \mathrm{for}\; s = 0, m \ge 0 \\
376  * {}_s Y_{\ell m}, & \mathrm{for}\; s > 0, m \ge 0 \\
377  * (-1)^{s + m} {}_s Y_{\ell m}, & \mathrm{for}\; s < 0, m < 0 \\
378  * (-1)^{m} {}_s Y_{\ell m}, & \mathrm{for}\; s = 0, m < 0 \\
379  * (-1)^{m + 1} {}_s Y_{\ell m}, & \mathrm{for}\; s > 0, m < 0
380  * \end{cases} \\
381  * &\equiv {}_s S_{\ell m}^{\mathrm{real}} {}_s Y_{\ell m},
382  * \f}
383  *
384  * where the unadorned \f${}_s Y_{\ell m}\f$ on the right-hand-sides follow the
385  * (older) convention of \cite Goldberg1966uu. Note the phase in these
386  * expressions is not completely standardized, so should be checked carefully
387  * whenever coefficient data is directly manipulated.
388  *
389  * For reference, we can compute the relation between Goldberg spin-weighted
390  * moments \f${}_s f_{\ell m}\f$, defined as
391  *
392  * \f[ {}_s f(\theta, \phi) = \sum_{\ell = 0}^{\ell_\mathrm{max}} \sum_{m =
393  * -\ell}^{\ell} {}_s f_{\ell m} {}_s Y_{\ell m}
394  * \f]
395  *
396  * so,
397  * \f[
398  * {}_s f_{\ell m} =
399  * \begin{cases}
400  * a_{\ell m}^{\mathrm{sharp}, \mathrm{real}} {}_s S_{\ell m}^{\mathrm{real}} +
401  * a_{\ell m}^{\mathrm{sharp}, \mathrm{imag}} {}_s S_{\ell m}^{\mathrm{imag}} &
402  * \mathrm{for} \; m \ge 0 \\
403  * \left(a_{\ell -m}^{\mathrm{sharp}, \mathrm{real}}\right)^* {}_s S_{\ell
404  * m}^{\mathrm{real}} + \left(a_{\ell -m}^{\mathrm{sharp},
405  * \mathrm{imag}}\right)^* {}_s S_{\ell m}^{\mathrm{imag}} &
406  * \mathrm{for} \; m < 0 \\
407  * \end{cases} \f]
408  *
409  *
410  * The resulting coefficients \f$a_{\ell m}\f$ are stored in a triangular,
411  * \f$\ell\f$-varies-fastest configuration. So, for instance, the first
412  * \f$\ell_\mathrm{max}\f$ entries contain the coefficients for \f$m=0\f$ and
413  * all \f$\ell\f$s, and the next \f$\ell_\mathrm{max} - 1\f$ entries contain the
414  * coefficients for \f$m=1\f$ and \f$1 \le \ell \le \ell_\mathrm{max} \f$, and
415  * so on.
416  */
417 template <typename TagList, ComplexRepresentation Representation =
418  ComplexRepresentation::Interleaved>
420 
421 /// \cond
422 template <typename... TransformTags, ComplexRepresentation Representation>
423 struct SwshTransform<tmpl::list<TransformTags...>, Representation> {
424  static constexpr int spin = db::const_item_type<
425  tmpl::front<tmpl::list<TransformTags...>>>::type::spin;
426 
427  static_assert(tmpl2::flat_all_v<
429  "All Tags in TagList submitted to SwshTransform must have the "
430  "same spin weight.");
431 
432  using return_tags = tmpl::list<Tags::SwshTransform<TransformTags>...>;
433  using argument_tags =
434  tmpl::list<TransformTags..., Tags::LMax, Tags::NumberOfRadialPoints>;
435 
436  static void apply(
437  const gsl::not_null<
439  const db::const_item_type<TransformTags>&... collocations,
440  const size_t l_max, const size_t number_of_radial_points) noexcept {
441  // forward to the version which takes parameter packs of vectors
442  apply_to_vectors(make_not_null(&get(*coefficients))...,
443  get(collocations)..., l_max, number_of_radial_points);
444  }
445 
446  // the convenience transform functions call the private member directly, so
447  // need to be friends
448  // redundant declaration due to forward declaration needs in earlier part of
449  // file and friend requirements
450  template <ComplexRepresentation FriendRepresentation, int FriendSpin,
451  typename... CoefficientThenCollocationTypes, size_t... Is>
452  // NOLINTNEXTLINE(readability-redundant-declaratioon)
453  friend void detail::swsh_transform_impl(
454  size_t, size_t, std::index_sequence<Is...>,
456  const CoefficientThenCollocationTypes&...) noexcept;
457 
458  // friend declaration for helper function from `SwshDerivatives.hpp`
459  template <typename Transform, typename FullTagList>
460  friend struct detail::dispatch_to_transform;
461 
462  private:
463  static void apply_to_vectors(
464  const gsl::not_null<typename db::item_type<
465  Tags::SwshTransform<TransformTags>>::type*>... coefficients,
466  const typename db::const_item_type<TransformTags>::type&... collocations,
467  size_t l_max, size_t number_of_radial_points) noexcept;
468 };
469 /// \endcond
470 
471 /*!
472  * \ingroup SwshGroup
473  * \brief A \ref DataBoxGroup mutate-compatible computational struct for
474  * performing several spin-weighted inverse spherical harmonic transforms.
475  * Internally dispatches to libsharp.
476  *
477  * \details
478  * Template Parameters:
479  * - `Representation`: The element of the `ComplexRepresentation` enum which
480  * parameterizes how the data passed to libsharp will be represented. Two
481  * options are available:
482  * - `ComplexRepresentation:Interleaved`: indicates that the real and imaginary
483  * parts of the collocation values will be passed to libsharp as pointers
484  * into existing complex data structures, minimizing copies, but maintaining
485  * a stride of 2 for 'adjacent' real or imaginary values.
486  * - `ComplexRepresentation::RealsThenImags`: indicates that the real and
487  * imaginary parts of the collocation values will be passed to libsharp as
488  * separate contiguous blocks. At current, this introduces both allocations
489  * and copies. **optimization note** In the future most of the allocations
490  * can be aggregated to calling code, which would pass in buffers by
491  * `not_null` pointers.
492  *
493  * For performance-sensitive code, both options should be tested, as each
494  * strategy has trade-offs.
495  * - `TagList`: A `tmpl::list` of Tags to be inverse transformed. The tags must
496  * represent the nodal data being transformed to.
497  *
498  * \see `SwshTransform` for mathematical notes regarding the libsharp modal
499  * representation taken as input for this computational struct.
500  */
501 template <typename TagList, ComplexRepresentation Representation =
502  ComplexRepresentation::Interleaved>
504 
505 /// \cond
506 template <typename... TransformTags, ComplexRepresentation Representation>
507 struct InverseSwshTransform<tmpl::list<TransformTags...>, Representation> {
508  static constexpr int spin = db::const_item_type<
509  tmpl::front<tmpl::list<TransformTags...>>>::type::spin;
510 
511  static_assert(
513  spin...>,
514  "All Tags in TagList submitted to InverseSwshTransform must have the "
515  "same spin weight.");
516 
517  using return_tags = tmpl::list<TransformTags...>;
518  using argument_tags = tmpl::list<Tags::SwshTransform<TransformTags>...,
519  Tags::LMax, Tags::NumberOfRadialPoints>;
520 
521  static void apply(
522  const gsl::not_null<db::item_type<TransformTags>*>... collocations,
523  const db::const_item_type<
524  Tags::SwshTransform<TransformTags>>&... coefficients,
525  const size_t l_max, const size_t number_of_radial_points) noexcept {
526  // forward to the version which takes parameter packs of vectors
527  apply_to_vectors(make_not_null(&get(*collocations))...,
528  get(coefficients)..., l_max, number_of_radial_points);
529  }
530 
531  // the convenience transform functions call the private member directly, so
532  // need to be friends
533  // redundant declaration due to forward declaration needs in earlier part of
534  // file and friend requirements
535  template <ComplexRepresentation FriendRepresentation, int FriendSpin,
536  typename... CollocationThenCoefficientTypes, size_t... Is>
537  // NOLINTNEXTLINE(readability-redundant-declaratioon)
538  friend void detail::inverse_swsh_transform_impl(
539  size_t, size_t, std::index_sequence<Is...>,
541  const CollocationThenCoefficientTypes&...) noexcept;
542 
543  // friend declaration for helper function from `SwshDerivatives.hpp`
544  template <typename Transform, typename FullTagList>
545  friend struct detail::dispatch_to_transform;
546 
547  private:
548  static void apply_to_vectors(
549  const gsl::not_null<
550  typename db::item_type<TransformTags>::type*>... collocations,
551  const typename db::const_item_type<
552  Tags::SwshTransform<TransformTags>>::type&... coefficients,
553  size_t l_max, size_t number_of_radial_points) noexcept;
554 };
555 
556 template <typename... TransformTags, ComplexRepresentation Representation>
557 void SwshTransform<tmpl::list<TransformTags...>, Representation>::
558  apply_to_vectors(
559  const gsl::not_null<typename db::item_type<
560  Tags::SwshTransform<TransformTags>>::type*>... coefficients,
561  const typename db::const_item_type<
562  TransformTags>::type&... collocations,
563  const size_t l_max, const size_t number_of_radial_points) noexcept {
564  EXPAND_PACK_LEFT_TO_RIGHT(coefficients->destructive_resize(
565  size_of_libsharp_coefficient_vector(l_max) * number_of_radial_points));
566 
567  // assemble a list of pointers into the collocation point data. This is
568  // required because libsharp expects pointers to pointers.
570  pre_transform_views.reserve(number_of_radial_points *
571  sizeof...(TransformTags));
572  std::vector<double*> pre_transform_collocation_data;
573  pre_transform_collocation_data.reserve(2 * number_of_radial_points *
574  sizeof...(TransformTags));
575 
576  // clang-tidy: const-cast, object is temporarily modified and returned to
577  // original state
578  EXPAND_PACK_LEFT_TO_RIGHT(detail::append_libsharp_collocation_pointers(
579  make_not_null(&pre_transform_collocation_data),
580  make_not_null(&pre_transform_views),
582  &const_cast<typename db::item_type<TransformTags>::type&>( // NOLINT
583  collocations)
584  .data()),
585  l_max, spin >= 0));
586 
587  std::vector<std::complex<double>*> post_transform_coefficient_data;
588  post_transform_coefficient_data.reserve(2 * number_of_radial_points *
589  sizeof...(TransformTags));
590 
591  EXPAND_PACK_LEFT_TO_RIGHT(detail::append_libsharp_coefficient_pointers(
592  make_not_null(&post_transform_coefficient_data),
593  make_not_null(&coefficients->data()), l_max));
594 
595  const size_t num_transforms =
596  (spin == 0 ? 2 : 1) * number_of_radial_points * sizeof...(TransformTags);
597  const auto* collocation_metadata =
598  &cached_collocation_metadata<Representation>(l_max);
599  const auto* alm_info =
600  cached_coefficients_metadata(l_max).get_sharp_alm_info();
601 
602  detail::execute_libsharp_transform_set(
603  SHARP_MAP2ALM, spin, make_not_null(&post_transform_coefficient_data),
604  make_not_null(&pre_transform_collocation_data),
605  make_not_null(collocation_metadata), alm_info, num_transforms);
606 
607  detail::conjugate_views<spin>(make_not_null(&pre_transform_views));
608 }
609 
610 template <typename... TransformTags, ComplexRepresentation Representation>
611 void InverseSwshTransform<tmpl::list<TransformTags...>, Representation>::
612  apply_to_vectors(
613  const gsl::not_null<
614  typename db::item_type<TransformTags>::type*>... collocations,
615  const typename db::const_item_type<
616  Tags::SwshTransform<TransformTags>>::type&... coefficients,
617  const size_t l_max, const size_t number_of_radial_points) noexcept {
618  EXPAND_PACK_LEFT_TO_RIGHT(collocations->destructive_resize(
619  number_of_swsh_collocation_points(l_max) * number_of_radial_points));
620 
621  std::vector<std::complex<double>*> pre_transform_coefficient_data;
622  pre_transform_coefficient_data.reserve(2 * number_of_radial_points *
623  sizeof...(TransformTags));
624  // clang-tidy: const-cast, object is temporarily modified and returned to
625  // original state
626  EXPAND_PACK_LEFT_TO_RIGHT(detail::append_libsharp_coefficient_pointers(
627  make_not_null(&pre_transform_coefficient_data),
629  &const_cast<typename db::item_type< // NOLINT
630  Tags::SwshTransform<TransformTags>>::type&>(coefficients)
631  .data()),
632  l_max));
633 
635  post_transform_views.reserve(number_of_radial_points *
636  sizeof...(TransformTags));
637  std::vector<double*> post_transform_collocation_data;
638  post_transform_collocation_data.reserve(2 * number_of_radial_points *
639  sizeof...(TransformTags));
640 
641  EXPAND_PACK_LEFT_TO_RIGHT(detail::append_libsharp_collocation_pointers(
642  make_not_null(&post_transform_collocation_data),
643  make_not_null(&post_transform_views),
644  make_not_null(&collocations->data()), l_max, true));
645 
646  // libsharp considers two arrays per transform when spin is not zero.
647  const size_t num_transforms =
648  (spin == 0 ? 2 : 1) * number_of_radial_points * sizeof...(TransformTags);
649  const auto* collocation_metadata =
650  &cached_collocation_metadata<Representation>(l_max);
651  const auto* alm_info =
652  cached_coefficients_metadata(l_max).get_sharp_alm_info();
653 
654  detail::execute_libsharp_transform_set(
655  SHARP_ALM2MAP, spin, make_not_null(&pre_transform_coefficient_data),
656  make_not_null(&post_transform_collocation_data),
657  make_not_null(collocation_metadata), alm_info, num_transforms);
658 
659  detail::conjugate_views<spin>(make_not_null(&post_transform_views));
660 
661  // The inverse transformed collocation data has just been placed in the
662  // memory blocks controlled by the `ComplexDataView`s. Finally, that data
663  // must be flushed back to the Variables.
664  for (auto& view : post_transform_views) {
665  view.copy_back_to_source();
666  }
667 }
668 /// \endcond
669 
670 namespace detail {
671 
672 template <ComplexRepresentation Representation, int Spin,
673  typename... ModalThenNodalTypes, size_t... Is>
674 void swsh_transform_impl(
675  const size_t l_max, const size_t number_of_radial_points,
678  first_coefficient,
679  const ModalThenNodalTypes&... coefficients_then_collocations) noexcept {
681  tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<Is, ComplexDataVector>,
683  apply_to_vectors(first_coefficient, coefficients_then_collocations...,
684  l_max, number_of_radial_points);
685 }
686 
687 template <ComplexRepresentation Representation, int Spin,
688  typename... NodalThenModalTypes, size_t... Is>
689 void inverse_swsh_transform_impl(
690  const size_t l_max, const size_t number_of_radial_points,
693  first_collocation,
694  const NodalThenModalTypes&... collocations_then_coefficients) noexcept {
696  tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<Is, ComplexDataVector>,
698  apply_to_vectors(first_collocation, collocations_then_coefficients...,
699  l_max, number_of_radial_points);
700 }
701 
702 // A metafunction for binning a provided tag list into `SwshTransform` objects
703 // according to spin-weight
704 template <int MinSpin, ComplexRepresentation Representation, typename TagList,
705  typename IndexSequence>
706 struct make_transform_list_impl;
707 
708 template <int MinSpin, ComplexRepresentation Representation, typename TagList,
709  int... Is>
710 struct make_transform_list_impl<MinSpin, Representation, TagList,
711  std::integer_sequence<int, Is...>> {
712  using type = tmpl::flatten<tmpl::list<tmpl::conditional_t<
713  not cpp17::is_same_v<get_tags_with_spin<Is + MinSpin, TagList>,
714  tmpl::list<>>,
716  tmpl::list<>>...>>;
717 };
718 
719 // A metafunction for binning a provided tag list into `InverseSwshTransform`
720 // objects according to spin-weight
721 template <int MinSpin, ComplexRepresentation Representation, typename TagList,
722  typename IndexSequence>
723 struct make_inverse_transform_list_impl;
724 
725 template <int MinSpin, ComplexRepresentation Representation, typename TagList,
726  int... Is>
727 struct make_inverse_transform_list_impl<MinSpin, Representation, TagList,
728  std::integer_sequence<int, Is...>> {
729  using type = tmpl::flatten<tmpl::list<tmpl::conditional_t<
730  not cpp17::is_same_v<get_tags_with_spin<Is + MinSpin, TagList>,
731  tmpl::list<>>,
733  Representation>,
734  tmpl::list<>>...>>;
735 };
736 } // namespace detail
737 
738 // @{
739 /// \ingroup SwshGroup
740 /// \brief Assemble a `tmpl::list` of `SwshTransform`s or
741 /// `InverseSwshTransform`s given a list of tags `TagList` that need to be
742 /// transformed. The `Representation` is the
743 /// `Spectral::Swsh::ComplexRepresentation` to use for the transformations.
744 ///
745 /// \details Up to five `SwshTransform`s or `InverseSwshTransform`s will be
746 /// returned, corresponding to the possible spin values. Any number of
747 /// transformations are aggregated into that set of `SwshTransform`s (or
748 /// `InverseSwshTransform`s). The number of transforms is up to five because the
749 /// libsharp utility only has capability to perform spin-weighted spherical
750 /// harmonic transforms for integer spin-weights from -2 to 2.
751 ///
752 /// \snippet Test_SwshTransform.cpp make_transform_list
753 template <ComplexRepresentation Representation, typename TagList>
754 using make_transform_list = typename detail::make_transform_list_impl<
755  -2, Representation, TagList,
757 
758 template <ComplexRepresentation Representation, typename TagList>
760  typename detail::make_inverse_transform_list_impl<
761  -2, Representation, TagList,
763 // @}
764 
765 /// \ingroup SwshGroup
766 /// \brief Assemble a `tmpl::list` of `SwshTransform`s given a list of
767 /// `Spectral::Swsh::Tags::Derivative<Tag, Derivative>` that need to be
768 /// computed. The `SwshTransform`s constructed by this type alias correspond to
769 /// the `Tag`s in the list.
770 ///
771 /// \details This is intended as a convenience utility for the first step of a
772 /// derivative routine, where one transforms the set of fields for which
773 /// derivatives are required.
774 ///
775 /// \snippet Test_SwshTransform.cpp make_transform_from_derivative_tags
776 template <ComplexRepresentation Representation, typename DerivativeTagList>
778  typename detail::make_transform_list_impl<
779  -2, Representation,
780  tmpl::transform<DerivativeTagList,
781  tmpl::bind<db::remove_tag_prefix, tmpl::_1>>,
783 } // namespace Swsh
784 } // namespace Spectral
constexpr size_t size_of_libsharp_coefficient_vector(const size_t l_max) noexcept
Convenience function for determining the number of spin-weighted spherical harmonics coefficients tha...
Definition: SwshCoefficients.hpp:34
A DataBox mutate-compatible computational struct for performing several spin-weighted inverse spheric...
Definition: SwshTransform.hpp:503
A DataBox mutate-compatible computational struct for performing several spin-weighted spherical harmo...
Definition: SwshTransform.hpp:419
typename detail::make_inverse_transform_list_impl< -2, Representation, TagList, decltype(std::make_integer_sequence< int, 5 >{})>::type make_inverse_transform_list
Assemble a tmpl::list of SwshTransforms or InverseSwshTransforms given a list of tags TagList that ne...
Definition: SwshTransform.hpp:762
#define EXPAND_PACK_LEFT_TO_RIGHT(...)
Expand a parameter pack evaluating the terms from left to right.
Definition: TMPL.hpp:562
constexpr size_t number_of_swsh_collocation_points(const size_t l_max) noexcept
Convenience function for determining the number of spin-weighted spherical harmonic collocation value...
Definition: SwshCollocation.hpp:25
constexpr bool flat_all_v
A non-short-circuiting logical AND between bools &#39;B"".
Definition: TMPL.hpp:505
Definition: Digraph.hpp:11
const CoefficientsMetadata & cached_coefficients_metadata(const size_t l_max) noexcept
Generation function for obtaining a CoefficientsMetadata object which is computed by the libsharp cal...
Definition: SwshCoefficients.cpp:32
Prefix tag representing the spin-weighted spherical harmonic transform of a spin-weighted scalar...
Definition: SwshTags.hpp:161
typename detail::make_transform_list_impl< -2, Representation, tmpl::transform< DerivativeTagList, tmpl::bind< db::remove_tag_prefix, tmpl::_1 > >, decltype(std::make_integer_sequence< int, 5 >{})>::type make_transform_list_from_derivative_tags
Assemble a tmpl::list of SwshTransforms given a list of Spectral::Swsh::Tags::Derivative<Tag, Derivative> that need to be computed. The SwshTransforms constructed by this type alias correspond to the Tags in the list.
Definition: SwshTransform.hpp:782
SpinWeighted< ComplexDataVector, Spin > inverse_swsh_transform(const size_t l_max, const size_t number_of_radial_points, const SpinWeighted< ComplexModalVector, Spin > &libsharp_coefficients) noexcept
Perform an inverse libsharp spin-weighted spherical harmonic transform on a single supplied SpinWeigh...
Definition: SwshTransform.cpp:118
SpinWeighted< ComplexModalVector, Spin > swsh_transform(const size_t l_max, const size_t number_of_radial_points, const SpinWeighted< ComplexDataVector, Spin > &collocation) noexcept
Perform a forward libsharp spin-weighted spherical harmonic transform on a single supplied SpinWeight...
Definition: SwshTransform.cpp:107
Definition: Determinant.hpp:11
Tag for the number of radial grid points in the three-dimensional representation of radially concentr...
Definition: SwshTags.hpp:186
Make a spin-weighted type T with spin-weight Spin. Mathematical operators are restricted to addition...
Definition: SpinWeighted.hpp:25
#define SPECTRE_ALWAYS_INLINE
Always inline a function. Only use this if you benchmarked the code.
Definition: ForceInline.hpp:16
Stores a collection of complex function values.
Definition: ComplexDataVector.hpp:47
A class for storing complex spectral coefficients on a spectral grid.
Definition: ComplexModalVector.hpp:39
Defines class Variables.
typename detail::make_transform_list_impl< -2, Representation, TagList, decltype(std::make_integer_sequence< int, 5 >{})>::type make_transform_list
Assemble a tmpl::list of SwshTransforms or InverseSwshTransforms given a list of tags TagList that ne...
Definition: SwshTransform.hpp:756
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args) noexcept
Apply the invokable f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1623
Tag for the maximum spin-weighted spherical harmonic l; sets angular resolution.
Definition: SwshTags.hpp:178
Wraps the template metaprogramming library used (brigand)
typename DataBox_detail::item_type_impl< TagList, Tag >::type item_type
Get the type that can be written to the Tag. If it is a base tag then a TagList must be passed as a s...
Definition: DataBoxTag.hpp:461
Defines functions and classes from the GSL.
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion, but it may be necessary to perform the conversion explicitly when type deduction is desired.
Definition: Gsl.hpp:879
Functionality associated with a particular choice of basis functions and quadrature for spectral oper...
Definition: Chebyshev.cpp:19
ComplexRepresentation
A set of labels for the possible representations of complex numbers for the ComplexDataView and the c...
Definition: ComplexDataView.hpp:57
Defines classes SimpleTag, PrefixTag, ComputeTag and several functions for retrieving tag info...
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12