SwshTestHelpers.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cmath>
7 #include <complex>
8 #include <cstddef>
9 #include <sharp_cxx.h>
10 
11 #include "DataStructures/ComplexModalVector.hpp"
12 #include "DataStructures/DataVector.hpp"
13 #include "NumericalAlgorithms/Spectral/ComplexDataView.hpp"
14 #include "NumericalAlgorithms/Spectral/SwshCoefficients.hpp"
15 #include "NumericalAlgorithms/Spectral/SwshCollocation.hpp"
16 #include "NumericalAlgorithms/Spectral/SwshTags.hpp" // IWYU pragma: keep
17 #include "Utilities/Gsl.hpp"
19 
20 /// \cond
21 class ComplexDataVector;
22 /// \endcond
23 
24 namespace Spectral {
25 namespace Swsh {
26 namespace TestHelpers {
27 
28 // returns the factorial of the argument as a double so that an approximate
29 // value can be given for larger input quantities. Note that the spin-weighted
30 // harmonic function requires a factorial of l + m, so harmonics above l~12
31 // would be untestable if the factorial only returned size_t's.
32 double factorial(size_t arg) noexcept;
33 
34 // Note that the methods for computing the spin-weighted spherical harmonics and
35 // their derivatives below are both 1) poorly optimized (they use many
36 // computations per grid point evaluated) and 2) not terribly accurate (the
37 // analytic expressions require evaluation of ratios of factorials, losing
38 // numerical precision rapidly). However, they are comparatively easy to
39 // manually check for correctness, which is critical to offer a reliable measure
40 // of the spin-weighted transforms.
41 
42 // Analytic form for the spin-weighted spherical harmonic function, for testing
43 // purposes. The formula is from [Goldberg
44 // et. al.](https://aip.scitation.org/doi/10.1063/1.1705135)
45 std::complex<double> spin_weighted_spherical_harmonic(int s, int l, int m,
46  double theta,
47  double phi) noexcept;
48 
49 // Returns the value of the spin-weighted derivative `DerivativeKind` of the
50 // spherical harmonic basis function \f${}_s Y_{l m}\f$ at angular location
51 // (`theta`, `phi`) using the recurrence identities for spin-weighted
52 // derivatives.
53 template <typename DerivativeKind>
54 std::complex<double> derivative_of_spin_weighted_spherical_harmonic(
55  int s, int l, int m, double theta, double phi) noexcept;
56 
57 // Performs the generation of spin-weighted coefficients into a supplied
58 // ComplexModalVector `to_fill`, passed by pointer. The resulting random
59 // coefficients must then be adjusted in two ways:
60 // - The modes l < |s| are always zero, so must have zero coefficient to be
61 // compatible with tests
62 // - Modes with m=0 obey reality conditions, so must also be adjusted
63 // The libsharp coefficient representation is not currently presented as an
64 // interface in SpECTRE. However, to help maintain the Swsh library, further
65 // notes on the coefficient representation can be found in the comments for the
66 // `Coefficients` class in `SwshCoefficients.hpp`.
67 template <int Spin, typename Distribution, typename Generator>
68 void generate_swsh_modes(const gsl::not_null<ComplexModalVector*> to_fill,
69  const gsl::not_null<Generator*> generator,
70  const gsl::not_null<Distribution*> distribution,
71  const size_t number_of_radial_points,
72  const size_t l_max) noexcept {
73  fill_with_random_values(to_fill, generator, distribution);
74 
75  auto spherical_harmonic_lm =
76  cached_coefficients_metadata(l_max).get_sharp_alm_info();
77  for (size_t i = 0; i < number_of_radial_points; i++) {
78  // adjust the m = 0 modes for the reality conditions in libsharp
79  for (size_t l = 0; l < l_max + 1; l++) {
80  (*to_fill)[l + i * size_of_libsharp_coefficient_vector(l_max)] = real(
81  (*to_fill)[l + i * size_of_libsharp_coefficient_vector(l_max) / 2]);
82  (*to_fill)[l +
83  (2 * i + 1) * size_of_libsharp_coefficient_vector(l_max) / 2] =
84  imag((*to_fill)[l + (2 * i + 1) *
86  2]);
87  }
88  for (size_t m = 0; m < static_cast<size_t>(spherical_harmonic_lm->nm);
89  ++m) {
90  for (size_t l = m; l <= l_max; ++l) {
91  // clang-tidy do not use pointer arithmetic
92  // pointer arithmetic here is unavoidable as we are interacting with
93  // the libsharp structures
94  const auto m_start =
95  static_cast<size_t>(spherical_harmonic_lm->mvstart[m]); // NOLINT
96  const auto l_stride =
97  static_cast<size_t>(spherical_harmonic_lm->stride);
98  // modes for l < |s| are zero
99  if (static_cast<int>(l) < abs(Spin)) {
100  (*to_fill)[m_start + l * l_stride +
101  i * size_of_libsharp_coefficient_vector(l_max)] = 0.0;
102  (*to_fill)[m_start + l * l_stride +
103  (2 * i + 1) * size_of_libsharp_coefficient_vector(l_max) /
104  2] = 0.0;
105  }
106  }
107  }
108  }
109 }
110 
111 // Computes a set of collocation points from a set of spin-weighted spherical
112 // harmonic coefficients. This function takes care of the nastiness associated
113 // with appropriately iterating over the sharp representation of coefficients
114 // and storing the result computed from the input `basis_function` evaluated
115 // with arguments (size_t s, size_t l, size_t m, double theta, double phi),
116 // which is compatible with the above `spin_weighted_spherical_harmonic`
117 // and `derivative_of_spin_weighted_spherical_harmonic` functions.
118 // This function does not have easy test verification, but acts as an
119 // independent computation of the spherical harmonic collocation values, so the
120 // agreement (from `Test_SwshTransformJob.cpp`) with the libsharp transform
121 // lends credibility to both methods.
122 template <int Spin, ComplexRepresentation Representation,
123  typename BasisFunction>
124 void swsh_collocation_from_coefficients_and_basis_func(
125  const gsl::not_null<ComplexDataVector*> collocation_data,
126  const ComplexModalVector& coefficient_data, const size_t l_max,
127  const size_t number_of_radial_points,
128  const BasisFunction basis_function) noexcept {
129  auto& spherical_harmonic_collocation =
130  cached_collocation_metadata<Representation>(l_max);
131  auto spherical_harmonic_lm =
132  cached_coefficients_metadata(l_max).get_sharp_alm_info();
133 
134  *collocation_data = 0.0;
135  for (size_t i = 0; i < number_of_radial_points; ++i) {
136  for (auto j : spherical_harmonic_collocation) {
137  for (size_t m = 0; m < static_cast<size_t>(spherical_harmonic_lm->nm);
138  ++m) {
139  for (size_t l = m; l <= l_max; ++l) {
140  // clang-tidy do not use pointer arithmetic
141  // pointer arithmetic here is unavoidable as we are interacting with
142  // the libsharp structures
143  const auto m_start =
144  static_cast<size_t>(spherical_harmonic_lm->mvstart[m]); // NOLINT
145  const auto l_stride =
146  static_cast<size_t>(spherical_harmonic_lm->stride);
147  // see the documentation associated with Swsh::TransformJob in
148  // SwshTransformJob.hpp for a detailed explanation of the libsharp
149  // collocation representation, and why there must be four steps with
150  // sign transformations.
151  (*collocation_data)[j.offset +
152  i * spherical_harmonic_collocation.size()] +=
153  sharp_swsh_sign(Spin, m, true) *
154  coefficient_data[m_start + l * l_stride +
156  basis_function(Spin, l, m, j.theta, j.phi);
157 
158  if (m != 0) {
159  (*collocation_data)[j.offset +
160  i * spherical_harmonic_collocation.size()] +=
161  sharp_swsh_sign(Spin, -m, true) *
162  conj(coefficient_data[m_start + l * l_stride +
164  l_max)]) *
165  basis_function(Spin, l, -m, j.theta, j.phi);
166  }
167  (*collocation_data)[j.offset +
168  i * spherical_harmonic_collocation.size()] +=
169  sharp_swsh_sign(Spin, m, false) * std::complex<double>(0.0, 1.0) *
170  coefficient_data[m_start + l * l_stride +
171  (2 * i + 1) *
173  2] *
174  basis_function(Spin, l, m, j.theta, j.phi);
175  if (m != 0) {
176  (*collocation_data)[j.offset +
177  i * spherical_harmonic_collocation.size()] +=
178  sharp_swsh_sign(Spin, -m, false) *
179  std::complex<double>(0.0, 1.0) *
180  conj(coefficient_data[m_start + l * l_stride +
181  (2 * i + 1) *
183  l_max) /
184  2]) *
185  basis_function(Spin, l, -m, j.theta, j.phi);
186  }
187  }
188  }
189  }
190  }
191 }
192 } // namespace TestHelpers
193 } // namespace Swsh
194 } // namespace Spectral
constexpr size_t size_of_libsharp_coefficient_vector(const size_t l_max) noexcept
Convenience function for determining the number of spin-weighted spherical harmonics coefficients tha...
Definition: SwshCoefficients.hpp:34
Definition: TestCreation.hpp:14
void phi(gsl::not_null< tnsr::iaa< DataType, SpatialDim, Frame > *> phi, const Scalar< DataType > &lapse, const tnsr::i< DataType, SpatialDim, Frame > &deriv_lapse, const tnsr::I< DataType, SpatialDim, Frame > &shift, const tnsr::iJ< DataType, SpatialDim, Frame > &deriv_shift, const tnsr::ii< DataType, SpatialDim, Frame > &spatial_metric, const tnsr::ijj< DataType, SpatialDim, Frame > &deriv_spatial_metric) noexcept
Computes the auxiliary variable used by the generalized harmonic formulation of Einstein&#39;s equations...
const CoefficientsMetadata & cached_coefficients_metadata(const size_t l_max) noexcept
Generation function for obtaining a CoefficientsMetadata object which is computed by the libsharp cal...
Definition: SwshCoefficients.cpp:32
Helper functions for data structures used in unit tests.
void fill_with_random_values(const gsl::not_null< T *> data, const gsl::not_null< UniformRandomBitGenerator *> generator, const gsl::not_null< RandomNumberDistribution *> distribution) noexcept
Fill an existing data structure with random values.
Definition: MakeWithRandomValues.hpp:144
Stores a collection of complex function values.
Definition: ComplexDataVector.hpp:47
A class for storing complex spectral coefficients on a spectral grid.
Definition: ComplexModalVector.hpp:39
constexpr uint64_t factorial(const uint64_t n) noexcept
Compute the factorial of .
Definition: ConstantExpressions.hpp:88
Defines functions and classes from the GSL.
constexpr double sharp_swsh_sign(const int spin_weight, const int m, const bool real) noexcept
Compute the sign change between the libsharp convention and the set of spin-weighted spherical harmon...
Definition: SwshCoefficients.hpp:113
Functionality associated with a particular choice of basis functions and quadrature for spectral oper...
Definition: Chebyshev.cpp:19
ComplexRepresentation
A set of labels for the possible representations of complex numbers for the ComplexDataView and the c...
Definition: ComplexDataView.hpp:57
Require a pointer to not be a nullptr
Definition: Gsl.hpp:182