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 =
427 : tmpl::list<TransformTags..., Tags::LMax, Tags::NumberOfRadialPoints>;
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 = tmpl::list<Tags::SwshTransform<TransformTags>...,
511 : Tags::LMax, Tags::NumberOfRadialPoints>;
512 :
513 : static void apply(
514 : const gsl::not_null<typename TransformTags::type*>... collocations,
515 : const typename Tags::SwshTransform<TransformTags>::type&... coefficients,
516 : const size_t l_max, const size_t number_of_radial_points) {
517 : // forward to the version which takes parameter packs of vectors
518 : apply_to_vectors(make_not_null(&get(*collocations))...,
519 : get(coefficients)..., l_max, number_of_radial_points);
520 : }
521 :
522 : // the convenience transform functions call the private member directly, so
523 : // need to be friends
524 : // redundant declaration due to forward declaration needs in earlier part of
525 : // file and friend requirements
526 : template <ComplexRepresentation FriendRepresentation, int FriendSpin,
527 : typename... CollocationThenCoefficientTypes, size_t... Is>
528 : // NOLINTNEXTLINE(readability-redundant-declaration)
529 : friend void detail::inverse_swsh_transform_impl(
530 : size_t, size_t, std::index_sequence<Is...>,
531 : gsl::not_null<SpinWeighted<ComplexDataVector, FriendSpin>*>,
532 : const CollocationThenCoefficientTypes&...);
533 :
534 : // friend declaration for helper function from `SwshDerivatives.hpp`
535 : template <typename Transform, typename FullTagList>
536 : friend struct detail::dispatch_to_transform;
537 :
538 : private:
539 : static void apply_to_vectors(
540 : const gsl::not_null<typename TransformTags::type::type*>... collocations,
541 : const typename Tags::SwshTransform<
542 : TransformTags>::type::type&... coefficients,
543 : size_t l_max, size_t number_of_radial_points);
544 : };
545 :
546 : template <typename... TransformTags, ComplexRepresentation Representation>
547 : void SwshTransform<tmpl::list<TransformTags...>, Representation>::
548 : apply_to_vectors(const gsl::not_null<typename Tags::SwshTransform<
549 : TransformTags>::type::type*>... coefficients,
550 : const typename TransformTags::type::type&... collocations,
551 : const size_t l_max, const size_t number_of_radial_points) {
552 : EXPAND_PACK_LEFT_TO_RIGHT(coefficients->destructive_resize(
553 : size_of_libsharp_coefficient_vector(l_max) * number_of_radial_points));
554 :
555 : // assemble a list of pointers into the collocation point data. This is
556 : // required because libsharp expects pointers to pointers.
557 : std::vector<detail::ComplexDataView<Representation>> pre_transform_views;
558 : pre_transform_views.reserve(number_of_radial_points *
559 : sizeof...(TransformTags));
560 : std::vector<double*> pre_transform_collocation_data;
561 : pre_transform_collocation_data.reserve(2 * number_of_radial_points *
562 : sizeof...(TransformTags));
563 :
564 : // clang-tidy: const-cast, object is temporarily modified and returned to
565 : // original state
566 : EXPAND_PACK_LEFT_TO_RIGHT(detail::append_libsharp_collocation_pointers(
567 : make_not_null(&pre_transform_collocation_data),
568 : make_not_null(&pre_transform_views),
569 : make_not_null(&const_cast<typename TransformTags::type::type&>( // NOLINT
570 : collocations)
571 : .data()),
572 : l_max, spin >= 0));
573 :
574 : std::vector<std::complex<double>*> post_transform_coefficient_data;
575 : post_transform_coefficient_data.reserve(2 * number_of_radial_points *
576 : sizeof...(TransformTags));
577 :
578 : EXPAND_PACK_LEFT_TO_RIGHT(detail::append_libsharp_coefficient_pointers(
579 : make_not_null(&post_transform_coefficient_data),
580 : make_not_null(&coefficients->data()), l_max));
581 :
582 : const size_t num_transforms =
583 : (spin == 0 ? 2 : 1) * number_of_radial_points * sizeof...(TransformTags);
584 : const auto* collocation_metadata =
585 : &cached_collocation_metadata<Representation>(l_max);
586 : const auto* alm_info =
587 : cached_coefficients_metadata(l_max).get_sharp_alm_info();
588 :
589 : detail::execute_libsharp_transform_set(
590 : SHARP_MAP2ALM, spin, make_not_null(&post_transform_coefficient_data),
591 : make_not_null(&pre_transform_collocation_data),
592 : make_not_null(collocation_metadata), alm_info, num_transforms);
593 :
594 : detail::conjugate_views<spin>(make_not_null(&pre_transform_views));
595 : }
596 :
597 : template <typename... TransformTags, ComplexRepresentation Representation>
598 : void InverseSwshTransform<tmpl::list<TransformTags...>, Representation>::
599 : apply_to_vectors(const gsl::not_null<
600 : typename TransformTags::type::type*>... collocations,
601 : const typename Tags::SwshTransform<
602 : TransformTags>::type::type&... coefficients,
603 : const size_t l_max, const size_t number_of_radial_points) {
604 : EXPAND_PACK_LEFT_TO_RIGHT(collocations->destructive_resize(
605 : number_of_swsh_collocation_points(l_max) * number_of_radial_points));
606 :
607 : std::vector<std::complex<double>*> pre_transform_coefficient_data;
608 : pre_transform_coefficient_data.reserve(2 * number_of_radial_points *
609 : sizeof...(TransformTags));
610 : // clang-tidy: const-cast, object is temporarily modified and returned to
611 : // original state
612 : EXPAND_PACK_LEFT_TO_RIGHT(detail::append_libsharp_coefficient_pointers(
613 : make_not_null(&pre_transform_coefficient_data),
614 : make_not_null(&const_cast<typename Tags::SwshTransform< // NOLINT
615 : TransformTags>::type::type&>(coefficients)
616 : .data()),
617 : l_max));
618 :
619 : std::vector<detail::ComplexDataView<Representation>> post_transform_views;
620 : post_transform_views.reserve(number_of_radial_points *
621 : sizeof...(TransformTags));
622 : std::vector<double*> post_transform_collocation_data;
623 : post_transform_collocation_data.reserve(2 * number_of_radial_points *
624 : sizeof...(TransformTags));
625 :
626 : EXPAND_PACK_LEFT_TO_RIGHT(detail::append_libsharp_collocation_pointers(
627 : make_not_null(&post_transform_collocation_data),
628 : make_not_null(&post_transform_views),
629 : make_not_null(&collocations->data()), l_max, true));
630 :
631 : // libsharp considers two arrays per transform when spin is not zero.
632 : const size_t num_transforms =
633 : (spin == 0 ? 2 : 1) * number_of_radial_points * sizeof...(TransformTags);
634 : const auto* collocation_metadata =
635 : &cached_collocation_metadata<Representation>(l_max);
636 : const auto* alm_info =
637 : cached_coefficients_metadata(l_max).get_sharp_alm_info();
638 :
639 : detail::execute_libsharp_transform_set(
640 : SHARP_ALM2MAP, spin, make_not_null(&pre_transform_coefficient_data),
641 : make_not_null(&post_transform_collocation_data),
642 : make_not_null(collocation_metadata), alm_info, num_transforms);
643 :
644 : detail::conjugate_views<spin>(make_not_null(&post_transform_views));
645 :
646 : // The inverse transformed collocation data has just been placed in the
647 : // memory blocks controlled by the `ComplexDataView`s. Finally, that data
648 : // must be flushed back to the Variables.
649 : for (auto& view : post_transform_views) {
650 : view.copy_back_to_source();
651 : }
652 : }
653 : /// \endcond
654 :
655 : namespace detail {
656 :
657 : template <ComplexRepresentation Representation, int Spin,
658 : typename... ModalThenNodalTypes, size_t... Is>
659 : void swsh_transform_impl(
660 : const size_t l_max, const size_t number_of_radial_points,
661 : std::index_sequence<Is...> /*meta*/,
662 : const gsl::not_null<SpinWeighted<ComplexModalVector, Spin>*>
663 : first_coefficient,
664 : const ModalThenNodalTypes&... coefficients_then_collocations) {
665 : SwshTransform<
666 : tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<Is, ComplexDataVector>,
667 : std::integral_constant<int, Spin>>...>>::
668 : apply_to_vectors(first_coefficient, coefficients_then_collocations...,
669 : l_max, number_of_radial_points);
670 : }
671 :
672 : template <ComplexRepresentation Representation, int Spin,
673 : typename... NodalThenModalTypes, size_t... Is>
674 : void inverse_swsh_transform_impl(
675 : const size_t l_max, const size_t number_of_radial_points,
676 : std::index_sequence<Is...> /*meta*/,
677 : const gsl::not_null<SpinWeighted<ComplexDataVector, Spin>*>
678 : first_collocation,
679 : const NodalThenModalTypes&... collocations_then_coefficients) {
680 : InverseSwshTransform<
681 : tmpl::list<::Tags::SpinWeighted<::Tags::TempScalar<Is, ComplexDataVector>,
682 : std::integral_constant<int, Spin>>...>>::
683 : apply_to_vectors(first_collocation, collocations_then_coefficients...,
684 : l_max, number_of_radial_points);
685 : }
686 :
687 : template <ComplexRepresentation Representation, typename PartitionedTagList>
688 : struct make_transform_list_impl;
689 :
690 : template <ComplexRepresentation Representation, typename... Partition>
691 : struct make_transform_list_impl<Representation, tmpl::list<Partition...>> {
692 : using type = tmpl::list<SwshTransform<Partition, Representation>...>;
693 : };
694 :
695 : template <ComplexRepresentation Representation, typename PartitionedTagList>
696 : struct make_inverse_transform_list_impl;
697 :
698 : template <ComplexRepresentation Representation, typename... Partition>
699 : struct make_inverse_transform_list_impl<Representation,
700 : tmpl::list<Partition...>> {
701 : using type = tmpl::list<InverseSwshTransform<Partition, Representation>...>;
702 : };
703 : } // namespace detail
704 :
705 : /// @{
706 : /// \ingroup SwshGroup
707 : /// \brief Assemble a `tmpl::list` of `SwshTransform`s or
708 : /// `InverseSwshTransform`s given a list of tags `TagList` that need to be
709 : /// transformed. The `Representation` is the
710 : /// `Spectral::Swsh::ComplexRepresentation` to use for the transformations.
711 : ///
712 : /// \snippet Test_SwshTransform.cpp make_transform_list
713 : template <ComplexRepresentation Representation, typename TagList>
714 1 : using make_transform_list = typename detail::make_transform_list_impl<
715 : Representation, partition_tags_by_spin<TagList>>::type;
716 :
717 : template <ComplexRepresentation Representation, typename TagList>
718 1 : using make_inverse_transform_list =
719 : typename detail::make_inverse_transform_list_impl<
720 : Representation, partition_tags_by_spin<TagList>>::type;
721 : /// @}
722 :
723 : /// \ingroup SwshGroup
724 : /// \brief Assemble a `tmpl::list` of `SwshTransform`s given a list of
725 : /// `Spectral::Swsh::Tags::Derivative<Tag, Derivative>` that need to be
726 : /// computed. The `SwshTransform`s constructed by this type alias correspond to
727 : /// the `Tag`s in the list.
728 : ///
729 : /// \details This is intended as a convenience utility for the first step of a
730 : /// derivative routine, where one transforms the set of fields for which
731 : /// derivatives are required.
732 : ///
733 : /// \snippet Test_SwshTransform.cpp make_transform_from_derivative_tags
734 : template <ComplexRepresentation Representation, typename DerivativeTagList>
735 1 : using make_transform_list_from_derivative_tags = make_transform_list<
736 : Representation,
737 : tmpl::remove_duplicates<tmpl::transform<
738 : DerivativeTagList, tmpl::bind<db::remove_tag_prefix, tmpl::_1>>>>;
739 :
740 : /// \ingroup SwshGroup
741 : /// \brief Convert spin-weighted spherical harmonic data to a new set of
742 : /// collocation points (either downsampling or upsampling)
743 : template <int Spin>
744 1 : void interpolate_to_collocation(
745 : gsl::not_null<SpinWeighted<ComplexDataVector, Spin>*> target,
746 : const SpinWeighted<ComplexDataVector, Spin>& source, size_t target_l_max,
747 : size_t source_l_max, size_t number_of_radial_points);
748 :
749 : } // namespace Swsh
750 : } // namespace Spectral
|
