YlmSpherepack.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <cmath>
8 #include <cstddef>
9 #include <utility>
10 #include <vector>
11 
12 #include "ApparentHorizons/YlmSpherepackHelper.hpp"
13 #include "DataStructures/DataVector.hpp"
15 #include "Utilities/Blas.hpp"
17 #include "Utilities/Gsl.hpp"
18 
19 /// \ingroup SpectralGroup
20 /// \brief C++ interface to SPHEREPACK.
22  public:
23  /// Type returned by gradient function.
24  using FirstDeriv = tnsr::i<DataVector, 2, Frame::Logical>;
25  /// Type returned by second derivative function.
26  using SecondDeriv = tnsr::ij<DataVector, 2, Frame::Logical>;
27 
28  /// Cached information to interpolate to the `target` point.
31  const std::array<double, 2>& target);
32  double cos_theta;
33  // cos(m*phi) and sin(m*phi)
34  std::vector<double> cos_m_phi, sin_m_phi;
35  // pbar_factor[m] = Pbar(m,m)*sin(theta)^m
36  std::vector<double> pbar_factor;
37  };
38 
39  /// Type to hold cached information at a set of target interpolation
40  /// points.
42 
43  /// Here l_max and m_max are the largest fully-represented l and m in
44  /// the Ylm expansion.
45  YlmSpherepack(size_t l_max, size_t m_max) noexcept;
46 
47  ///@{
48  /// Static functions to return the correct sizes of vectors of
49  /// collocation points and spectral coefficients for a given l_max
50  /// and m_max. Useful for allocating space without having to create
51  /// a YlmSpherepack.
52  SPECTRE_ALWAYS_INLINE static constexpr size_t physical_size(
53  const size_t l_max, const size_t m_max) noexcept {
54  return (l_max + 1) * (2 * m_max + 1);
55  }
56  /// \note `spectral_size` is the size of the buffer that holds the
57  /// coefficients; it is not the number of coefficients (which is
58  /// \f$m_{\rm max}^2+(l_{\rm max}-m_{\rm max})(2m_{\rm max}+1)\f$).
59  /// To simplify its internal indexing, SPHEREPACK uses a buffer with
60  /// more space than necessary. See SpherepackIterator for
61  /// how to index the coefficients in the buffer.
62  SPECTRE_ALWAYS_INLINE static constexpr size_t spectral_size(
63  const size_t l_max, const size_t m_max) noexcept {
64  return 2 * (l_max + 1) * (m_max + 1);
65  }
66  ///@}
67 
68  ///@{
69  /// Sizes in physical and spectral space for this instance.
70  size_t l_max() const noexcept { return l_max_; }
71  size_t m_max() const noexcept { return m_max_; }
72  size_t physical_size() const noexcept { return n_theta_ * n_phi_; }
73  size_t spectral_size() const noexcept { return spectral_size_; }
74  ///@}
75 
76  std::array<size_t, 2> physical_extents() const noexcept {
77  return {{n_theta_, n_phi_}};
78  }
79 
80  ///@{
81  /// Collocation points theta and phi.
82  const std::vector<double>& theta_points() const noexcept;
83  const std::vector<double>& phi_points() const noexcept;
85  ///@}
86 
87  ///@{
88  /// Spectral transformations.
89  /// To act on a slice of the input and output arrays, specify strides
90  /// and offsets.
91  void phys_to_spec(gsl::not_null<double*> spectral_coefs,
92  gsl::not_null<const double*> collocation_values,
93  size_t physical_stride = 1, size_t physical_offset = 0,
94  size_t spectral_stride = 1,
95  size_t spectral_offset = 0) const noexcept {
96  phys_to_spec_impl(spectral_coefs, collocation_values, physical_stride,
97  physical_offset, spectral_stride, spectral_offset, false);
98  }
99  void spec_to_phys(gsl::not_null<double*> collocation_values,
100  gsl::not_null<const double*> spectral_coefs,
101  size_t spectral_stride = 1, size_t spectral_offset = 0,
102  size_t physical_stride = 1,
103  size_t physical_offset = 0) const noexcept {
104  spec_to_phys_impl(collocation_values, spectral_coefs, spectral_stride,
105  spectral_offset, physical_stride, physical_offset, false);
106  };
107  ///@}
108 
109  ///@{
110  /// Spectral transformations where `collocation_values` and
111  /// `spectral_coefs` are assumed to point to 3-dimensional arrays
112  /// (I1 x S2 topology), and the transformations are done for all
113  /// 'radial' points at once by internally looping over all values of
114  /// the offset from zero to `stride`-1 (the physical and spectral
115  /// strides are equal and are called `stride`).
117  gsl::not_null<const double*> collocation_values,
118  size_t stride) const noexcept {
119  phys_to_spec_impl(spectral_coefs, collocation_values, stride, 0, stride, 0,
120  true);
121  }
123  gsl::not_null<const double*> spectral_coefs,
124  size_t stride) const noexcept {
125  spec_to_phys_impl(collocation_values, spectral_coefs, stride, 0, stride, 0,
126  true);
127  };
128  ///@}
129 
130  ///@{
131  /// Simpler, less general interfaces to `phys_to_spec` and `spec_to_phys`.
132  /// Acts on a slice of the input and returns a unit-stride result.
133  DataVector phys_to_spec(const DataVector& collocation_values,
134  size_t physical_stride = 1,
135  size_t physical_offset = 0) const noexcept;
136  DataVector spec_to_phys(const DataVector& spectral_coefs,
137  size_t spectral_stride = 1,
138  size_t spectral_offset = 0) const noexcept;
139  ///@}
140 
141  ///@{
142  /// Simpler, less general interfaces to `phys_to_spec_all_offsets`
143  /// and `spec_to_phys_all_offsets`. Result has the same stride as
144  /// the input.
145  DataVector phys_to_spec_all_offsets(const DataVector& collocation_values,
146  size_t stride) const noexcept;
147  DataVector spec_to_phys_all_offsets(const DataVector& spectral_coefs,
148  size_t stride) const noexcept;
149  ///@}
150 
151  /// Computes Pfaffian derivative (df/dtheta, csc(theta) df/dphi) at
152  /// the collocation values.
153  /// To act on a slice of the input and output arrays, specify stride
154  /// and offset (assumed to be the same for input and output).
155  void gradient(const std::array<double*, 2>& df,
156  gsl::not_null<const double*> collocation_values,
157  size_t physical_stride = 1, size_t physical_offset = 0) const
158  noexcept;
159 
160  /// Same as `gradient`, but takes the spectral coefficients (rather
161  /// than collocation values) of the function. This is more
162  /// efficient if one happens to already have the spectral
163  /// coefficients.
164  /// To act on a slice of the input and output arrays, specify strides
165  /// and offsets.
167  gsl::not_null<const double*> spectral_coefs,
168  size_t spectral_stride = 1,
169  size_t spectral_offset = 0,
170  size_t physical_stride = 1,
171  size_t physical_offset = 0) const noexcept {
172  gradient_from_coefs_impl(df, spectral_coefs, spectral_stride,
173  spectral_offset, physical_stride, physical_offset,
174  false);
175  }
176 
177  ///@{
178  /// Same as `gradient` but pointers are assumed to point to
179  /// 3-dimensional arrays (I1 x S2 topology), and the gradient is
180  /// done for all 'radial' points at once by internally looping
181  /// over all values of the offset from zero to `stride`-1.
183  gsl::not_null<const double*> collocation_values,
184  size_t stride = 1) const noexcept;
185 
187  const std::array<double*, 2>& df,
188  gsl::not_null<const double*> spectral_coefs, size_t stride = 1) const
189  noexcept {
190  gradient_from_coefs_impl(df, spectral_coefs, stride, 0, stride, 0, true);
191  }
192  ///@}
193 
194  ///@{
195  /// Simpler, less general interfaces to `gradient`.
196  /// Acts on a slice of the input and returns a unit-stride result.
197  FirstDeriv gradient(const DataVector& collocation_values,
198  size_t physical_stride = 1,
199  size_t physical_offset = 0) const noexcept;
200  FirstDeriv gradient_from_coefs(const DataVector& spectral_coefs,
201  size_t spectral_stride = 1,
202  size_t spectral_offset = 0) const noexcept;
203  ///@}
204 
205  ///@{
206  /// Simpler, less general interfaces to `gradient_all_offsets`.
207  /// Result has the same stride as the input.
208  FirstDeriv gradient_all_offsets(const DataVector& collocation_values,
209  size_t stride = 1) const noexcept;
211  size_t stride = 1) const noexcept;
212  ///@}
213 
214  /// Computes Laplacian in physical space.
215  /// To act on a slice of the input and output arrays, specify stride
216  /// and offset (assumed to be the same for input and output).
218  gsl::not_null<const double*> collocation_values,
219  size_t physical_stride = 1,
220  size_t physical_offset = 0) const noexcept;
221 
222  /// Same as `scalar_laplacian` above, but the input is the spectral
223  /// coefficients (rather than collocation values) of the function.
224  /// This is more efficient if one happens to already have the
225  /// spectral coefficients.
226  /// To act on a slice of the input and output arrays, specify strides
227  /// and offsets.
229  gsl::not_null<const double*> spectral_coefs,
230  size_t spectral_stride = 1,
231  size_t spectral_offset = 0,
232  size_t physical_stride = 1,
233  size_t physical_offset = 0) const noexcept;
234 
235  ///@{
236  /// Simpler, less general interfaces to `scalar_laplacian`.
237  /// Acts on a slice of the input and returns a unit-stride result.
238  DataVector scalar_laplacian(const DataVector& collocation_values,
239  size_t physical_stride = 1,
240  size_t physical_offset = 0) const noexcept;
241  DataVector scalar_laplacian_from_coefs(const DataVector& spectral_coefs,
242  size_t spectral_stride = 1,
243  size_t spectral_offset = 0) const
244  noexcept;
245  ///@}
246 
247  /// Computes Pfaffian first and second derivative in physical space.
248  /// The first derivative is \f$df(i) = d_i f\f$, and the
249  /// second derivative is \f$ddf(i,j) = d_i (d_j f)\f$,
250  /// where \f$d_0 = d/d\theta\f$ and \f$d_1 = csc(\theta) d/d\phi\f$.
251  /// ddf is not symmetric.
252  /// To act on a slice of the input and output arrays, specify stride
253  /// and offset (assumed to be the same for input and output).
256  gsl::not_null<const double*> collocation_values,
257  size_t physical_stride = 1,
258  size_t physical_offset = 0) const noexcept;
259 
260  /// Simpler, less general interface to second_derivative
262  const DataVector& collocation_values) const noexcept;
263 
264  /// Computes the integral over the sphere.
266  gsl::not_null<const double*> collocation_values,
267  size_t physical_stride = 1, size_t physical_offset = 0) const noexcept {
268  // clang-tidy: 'do not use pointer arithmetic'
269  return ddot_(n_theta_ * n_phi_, storage_.quadrature_weights.data(), 1,
270  collocation_values.get() + physical_offset, // NOLINT
271  physical_stride);
272  }
273 
274  /// Returns weights \f$w_i\f$ such that \f$sum_i (c_i w_i)\f$
275  /// is the definite integral, where \f$c_i\f$ are collocation values
276  /// at point i.
278  noexcept {
279  return storage_.quadrature_weights;
280  }
281 
282  /// Adds a constant (i.e. \f$f(\theta,\phi)\f$ += \f$c\f$) to the function
283  /// given by the spectral coefficients, by modifying the coefficients.
285  const gsl::not_null<DataVector*> spectral_coefs, const double c) const
286  noexcept {
287  // The factor of sqrt(8) is because of the normalization of
288  // SPHEREPACK's coefficients.
289  (*spectral_coefs)[0] += sqrt(8.0) * c;
290  }
291 
292  /// Returns the average of \f$f(\theta,\phi)\f$ over \f$(\theta,\phi)\f$.
293  SPECTRE_ALWAYS_INLINE double average(const DataVector& spectral_coefs) const
294  noexcept {
295  // The factor of sqrt(8) is because of the normalization of
296  // SPHEREPACK's coefficients. All other coefficients average to zero.
297  return spectral_coefs[0] / sqrt(8.0);
298  }
299 
300  /// Sets up the `InterpolationInfo` structure for interpolating onto
301  /// a set of target \f$(\theta,\phi)\f$ points. Does not depend on
302  /// the function being interpolated.
304  const std::vector<std::array<double, 2>>& target_points) const noexcept;
305 
306  /// Interpolates from `collocation_values` onto the points that have
307  /// been passed into the `set_up_interpolation_info` function.
308  /// To interpolate a different function on the same spectral grid, there
309  /// is no need to recompute `interpolation_info`.
310  /// If you specify stride and offset, acts on a slice of the input values.
311  /// The output has unit stride.
313  gsl::not_null<const double*> collocation_values,
314  const InterpolationInfo& interpolation_info,
315  size_t physical_stride = 1, size_t physical_offset = 0) const
316  noexcept;
317 
318  /// Same as `interpolate`, but assumes you have spectral coefficients.
319  /// This is more efficient if you already have the spectral coefficients
320  /// available.
321  /// If you specify stride and offset, acts on a slice of the input coefs.
322  /// The output has unit stride.
323  template <typename T>
325  const T& spectral_coefs,
326  const InterpolationInfo& interpolation_info,
327  size_t spectral_stride = 1,
328  size_t spectral_offset = 0) const noexcept;
329 
330  /// Simpler interface to `interpolate`. If you need to call this
331  /// repeatedly on different `spectral_coefs` or `collocation_values`
332  /// for the same target points, this is inefficient; instead use
333  /// `set_up_interpolation_info` and the functions that use
334  /// `InterpolationInfo`.
336  const DataVector& collocation_values,
337  const std::vector<std::array<double, 2>>& target_points) const noexcept;
339  const DataVector& spectral_coefs,
340  const std::vector<std::array<double, 2>>& target_points) const noexcept;
341 
342  /// Takes spectral coefficients compatible with `*this`, and either
343  /// prolongs them or restricts them to be compatible with `target`.
344  /// This is done by truncation (restriction) or padding with zeros
345  /// (prolongation).
346  DataVector prolong_or_restrict(const DataVector& spectral_coefs,
347  const YlmSpherepack& target) const noexcept;
348 
349  private:
350  // Spectral transformations and gradient.
351  // If `loop_over_offset` is true, then `collocation_values` and
352  // `spectral_coefs` are assumed to point to 3-dimensional
353  // arrays (I1 x S2 topology), and the transformations are done for
354  // all 'radial' points at once by looping over all values of the
355  // offset from zero to stride-1. If `loop_over_offset` is true,
356  // `physical_stride` must equal `spectral_stride`.
357  void phys_to_spec_impl(gsl::not_null<double*> spectral_coefs,
358  gsl::not_null<const double*> collocation_values,
359  size_t physical_stride = 1, size_t physical_offset = 0,
360  size_t spectral_stride = 1, size_t spectral_offset = 0,
361  bool loop_over_offset = false) const noexcept;
362  void spec_to_phys_impl(gsl::not_null<double*> collocation_values,
363  gsl::not_null<const double*> spectral_coefs,
364  size_t spectral_stride = 1, size_t spectral_offset = 0,
365  size_t physical_stride = 1, size_t physical_offset = 0,
366  bool loop_over_offset = false) const noexcept;
367  void gradient_from_coefs_impl(const std::array<double*, 2>& df,
368  gsl::not_null<const double*> spectral_coefs,
369  size_t spectral_stride = 1,
370  size_t spectral_offset = 0,
371  size_t physical_stride = 1,
372  size_t physical_offset = 0,
373  bool loop_over_offset = false) const noexcept;
374  void fill_scalar_work_arrays() const noexcept;
375  void fill_vector_work_arrays() const noexcept;
376  size_t l_max_, m_max_, n_theta_, n_phi_;
377  size_t spectral_size_;
378  mutable YlmSpherepack_detail::MemoryPool memory_pool_;
379  mutable YlmSpherepack_detail::Storage storage_;
380 }; // class YlmSpherepack
381 
382 bool operator==(const YlmSpherepack& lhs, const YlmSpherepack& rhs) noexcept;
383 bool operator!=(const YlmSpherepack& lhs, const YlmSpherepack& rhs) noexcept;
void gradient(const std::array< double *, 2 > &df, gsl::not_null< const double *> collocation_values, size_t physical_stride=1, size_t physical_offset=0) const noexcept
Computes Pfaffian derivative (df/dtheta, csc(theta) df/dphi) at the collocation values. To act on a slice of the input and output arrays, specify stride and offset (assumed to be the same for input and output).
const std::vector< double > & integration_weights() const noexcept
Returns weights such that is the definite integral, where are collocation values at point i...
Definition: YlmSpherepack.hpp:277
tnsr::i< DataVector, 2, Frame::Logical > FirstDeriv
Type returned by gradient function.
Definition: YlmSpherepack.hpp:24
void add_constant(const gsl::not_null< DataVector *> spectral_coefs, const double c) const noexcept
Adds a constant (i.e. += ) to the function given by the spectral coefficients, by modifying the coef...
Definition: YlmSpherepack.hpp:284
Cached information to interpolate to the target point.
Definition: YlmSpherepack.hpp:29
size_t l_max() const noexcept
Sizes in physical and spectral space for this instance.
Definition: YlmSpherepack.hpp:70
void scalar_laplacian_from_coefs(gsl::not_null< double *> scalar_laplacian, gsl::not_null< const double *> spectral_coefs, size_t spectral_stride=1, size_t spectral_offset=0, size_t physical_stride=1, size_t physical_offset=0) const noexcept
Same as scalar_laplacian above, but the input is the spectral coefficients (rather than collocation v...
void scalar_laplacian(gsl::not_null< double *> scalar_laplacian, gsl::not_null< const double *> collocation_values, size_t physical_stride=1, size_t physical_offset=0) const noexcept
Computes Laplacian in physical space. To act on a slice of the input and output arrays, specify stride and offset (assumed to be the same for input and output).
const std::vector< double > & theta_points() const noexcept
Collocation points theta and phi.
Definition: YlmSpherepack.cpp:391
double ddot_(const size_t &N, const double *X, const size_t &INCX, const double *Y, const size_t &INCY)
Definition: Blas.hpp:46
DataVector prolong_or_restrict(const DataVector &spectral_coefs, const YlmSpherepack &target) const noexcept
Takes spectral coefficients compatible with *this, and either prolongs them or restricts them to be c...
Definition: YlmSpherepack.cpp:890
void interpolate(gsl::not_null< std::vector< double > *> result, gsl::not_null< const double *> collocation_values, const InterpolationInfo &interpolation_info, size_t physical_stride=1, size_t physical_offset=0) const noexcept
Interpolates from collocation_values onto the points that have been passed into the set_up_interpolat...
double average(const DataVector &spectral_coefs) const noexcept
Returns the average of over .
Definition: YlmSpherepack.hpp:293
size_t m_max() const noexcept
Sizes in physical and spectral space for this instance.
Definition: YlmSpherepack.hpp:71
void phys_to_spec(gsl::not_null< double *> spectral_coefs, gsl::not_null< const double *> collocation_values, size_t physical_stride=1, size_t physical_offset=0, size_t spectral_stride=1, size_t spectral_offset=0) const noexcept
Spectral transformations. To act on a slice of the input and output arrays, specify strides and offse...
Definition: YlmSpherepack.hpp:91
YlmSpherepack(size_t l_max, size_t m_max) noexcept
Here l_max and m_max are the largest fully-represented l and m in the Ylm expansion.
Definition: YlmSpherepack.cpp:46
const std::vector< double > & phi_points() const noexcept
Collocation points theta and phi.
Definition: YlmSpherepack.cpp:409
C++ interface to SPHEREPACK.
Definition: YlmSpherepack.hpp:21
static constexpr size_t physical_size(const size_t l_max, const size_t m_max) noexcept
Static functions to return the correct sizes of vectors of collocation points and spectral coefficien...
Definition: YlmSpherepack.hpp:52
void phys_to_spec_all_offsets(gsl::not_null< double *> spectral_coefs, gsl::not_null< const double *> collocation_values, size_t stride) const noexcept
Spectral transformations where collocation_values and spectral_coefs are assumed to point to 3-dimens...
Definition: YlmSpherepack.hpp:116
void gradient_all_offsets(const std::array< double *, 2 > &df, gsl::not_null< const double *> collocation_values, size_t stride=1) const noexcept
Same as gradient but pointers are assumed to point to 3-dimensional arrays (I1 x S2 topology)...
double definite_integral(gsl::not_null< const double *> collocation_values, size_t physical_stride=1, size_t physical_offset=0) const noexcept
Computes the integral over the sphere.
Definition: YlmSpherepack.hpp:265
size_t spectral_size() const noexcept
Sizes in physical and spectral space for this instance.
Definition: YlmSpherepack.hpp:73
void gradient_from_coefs_all_offsets(const std::array< double *, 2 > &df, gsl::not_null< const double *> spectral_coefs, size_t stride=1) const noexcept
Same as gradient but pointers are assumed to point to 3-dimensional arrays (I1 x S2 topology)...
Definition: YlmSpherepack.hpp:186
#define SPECTRE_ALWAYS_INLINE
Always inline a function. Only use this if you benchmarked the code.
Definition: ForceInline.hpp:20
void second_derivative(const std::array< double *, 2 > &df, gsl::not_null< SecondDeriv *> ddf, gsl::not_null< const double *> collocation_values, size_t physical_stride=1, size_t physical_offset=0) const noexcept
Computes Pfaffian first and second derivative in physical space. The first derivative is ...
Definition: YlmSpherepack.cpp:423
tnsr::ij< DataVector, 2, Frame::Logical > SecondDeriv
Type returned by second derivative function.
Definition: YlmSpherepack.hpp:26
Defines macro to always inline a function.
Declares the interfaces for the BLAS used.
Defines a list of useful type aliases for tensors.
std::array< DataVector, 2 > theta_phi_points() const noexcept
Collocation points theta and phi.
Definition: YlmSpherepack.cpp:377
Stores a collection of function values.
Definition: DataVector.hpp:46
Defines functions and classes from the GSL.
InterpolationInfo set_up_interpolation_info(const std::vector< std::array< double, 2 >> &target_points) const noexcept
Sets up the InterpolationInfo structure for interpolating onto a set of target points. Does not depend on the function being interpolated.
Definition: YlmSpherepack.cpp:564
void interpolate_from_coefs(gsl::not_null< std::vector< double > *> result, const T &spectral_coefs, const InterpolationInfo &interpolation_info, size_t spectral_stride=1, size_t spectral_offset=0) const noexcept
Same as interpolate, but assumes you have spectral coefficients. This is more efficient if you alread...
void spec_to_phys(gsl::not_null< double *> collocation_values, gsl::not_null< const double *> spectral_coefs, size_t spectral_stride=1, size_t spectral_offset=0, size_t physical_stride=1, size_t physical_offset=0) const noexcept
Spectral transformations. To act on a slice of the input and output arrays, specify strides and offse...
Definition: YlmSpherepack.hpp:99
std::pair< FirstDeriv, SecondDeriv > first_and_second_derivative(const DataVector &collocation_values) const noexcept
Simpler, less general interface to second_derivative.
Definition: YlmSpherepack.cpp:514
void spec_to_phys_all_offsets(gsl::not_null< double *> collocation_values, gsl::not_null< const double *> spectral_coefs, size_t stride) const noexcept
Spectral transformations where collocation_values and spectral_coefs are assumed to point to 3-dimens...
Definition: YlmSpherepack.hpp:122
size_t physical_size() const noexcept
Sizes in physical and spectral space for this instance.
Definition: YlmSpherepack.hpp:72
static constexpr size_t spectral_size(const size_t l_max, const size_t m_max) noexcept
Definition: YlmSpherepack.hpp:62
void gradient_from_coefs(const std::array< double *, 2 > &df, gsl::not_null< const double *> spectral_coefs, size_t spectral_stride=1, size_t spectral_offset=0, size_t physical_stride=1, size_t physical_offset=0) const noexcept
Same as gradient, but takes the spectral coefficients (rather than collocation values) of the functio...
Definition: YlmSpherepack.hpp:166
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12