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