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