Line data Source code
1 0 : // 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 "DataStructures/DataVector.hpp"
13 : #include "DataStructures/DynamicBuffer.hpp"
14 : #include "DataStructures/Tensor/TypeAliases.hpp"
15 : #include "NumericalAlgorithms/SphericalHarmonics/SpherepackHelper.hpp"
16 : #include "Utilities/Blas.hpp"
17 : #include "Utilities/ForceInline.hpp"
18 : #include "Utilities/Gsl.hpp"
19 :
20 : /// Items related to spherical harmonics
21 : namespace ylm {
22 :
23 : /*!
24 : * \ingroup SpectralGroup
25 : *
26 : * \brief Defines the C++ interface to SPHEREPACK.
27 : *
28 : * \details The class `Spherepack` defines the C++ interface to the fortran
29 : * library SPHEREPACK used for computations on the surface of a sphere.
30 : *
31 : * Given a real-valued, scalar function \f$g(\theta, \phi)\f$, SPHEREPACK
32 : * expands it as:
33 : *
34 : * \f{align}
35 : * g(\theta, \phi)
36 : * &=\frac{1}{2}\sum_{l=0}^{l_{\max}}\bar P_l^0(\cos\theta) a_{l0}
37 : * +\sum_{l=1}^{l_{\max}}\sum_{m=1}^{\min(l, m_{\max})}\bar P_l^m(\cos\theta)\{
38 : * a_{lm}\cos m\phi -b_{lm}\sin m\phi\}\label{eq:spherepack_expansion}
39 : * \f}
40 : *
41 : * where \f$a_{lm}\f$ and \f$b_{lm}\f$ are real-valued
42 : * spectral coefficient arrays used by
43 : * SPHEREPACK, \f$P_l^m(x)\f$ are defined as
44 : *
45 : * \f{align}
46 : * \bar P_l^m(x)&=\sqrt{\frac{(2l+1)(l-m)!}{2(l+m)!}}\;P_{lm}(x)
47 : * \f}
48 : *
49 : * and \f$P_{nm}(x)\f$ are the associated Legendre polynomials as defined,
50 : * for example, in Jackson's "Classical Electrodynamics".
51 : *
52 : * #### Relationship to standard spherical harmonics
53 : *
54 : * The standard expansion of \f$g(\theta, \phi)\f$ in terms of scalar
55 : * spherical harmonics is
56 : * \f{align}
57 : * g(\theta, \phi)
58 : * &=
59 : * \sum_{l=0}^{l_{\max}}\sum_{m=-\min(l, m_{\max})}^{\min(l, m_{\max})}
60 : * A_{lm} Y_{lm}(\theta,\phi),
61 : * \f}
62 : * where \f$Y_{lm}(\theta,\phi)\f$ are the usual complex-valued scalar
63 : * spherical harmonics (as defined, for example, in
64 : * Jackson's "Classical Electrodynamics")
65 : * and \f$A_{lm}\f$ are complex coefficients.
66 : *
67 : * The relationship between the complex coefficients \f$A_{lm}\f$ and
68 : * SPHEREPACK's real-valued \f$a_{lm}\f$ and \f$b_{lm}\f$ is
69 : * \f{align}
70 : * a_{l0} & = \sqrt{\frac{2}{\pi}}A_{l0}&\qquad l\geq 0,\\
71 : * a_{lm} & = (-1)^m\sqrt{\frac{2}{\pi}} \mathrm{Re}(A_{lm})
72 : * &\qquad l\geq 1, m\geq 1, \\
73 : * b_{lm} & = (-1)^m\sqrt{\frac{2}{\pi}} \mathrm{Im}(A_{lm})
74 : * &\qquad l\geq 1, m\geq 1.
75 : * \f}
76 : *
77 : * \note If \f$g\f$ is real,
78 : * \f$A_{lm} = (-1)^m A^\star_{l -m}\f$ (where \f${}^\star\f$ means
79 : * a complex conjugate); this is why we don't need to consider \f$m<0\f$
80 : * in the previous formulas or in SPHEREPACK's expansion.
81 : *
82 : * #### Relationship to real-valued spherical harmonics
83 : *
84 : * Sometimes it is useful to expand a real-valued function in the form
85 : * \f{align}
86 : * g(\theta, \phi)
87 : * &= \sum_{l=0}^\infty\sum_{m=0}^l
88 : * \left[
89 : * c_{lm}\mathrm{Re}(Y_{lm}(\theta, \phi))+
90 : * d_{nm}\mathrm{Im}(Y_{lm}(\theta, \phi))
91 : * \right].
92 : * \f}
93 : * The coefficients here are therefore
94 : * \f{align}
95 : * c_{l0} &= A_{l0},\\
96 : * c_{lm} &= 2\mathrm{Re}(A_{lm}) \qquad m\geq 1,\\
97 : * d_{lm} &=-2\mathrm{Im}(A_{lm}).
98 : * \f}
99 : *
100 : * #### Modal and nodal representations
101 : *
102 : * Internally, SPHEREPACK can represent its expansion in two ways which we
103 : * will refer to as modal and nodal representations:
104 : *
105 : * -# modal: The spectral coefficient arrays \f$a_{lm}\f$ and \f$b_{lm}\f$,
106 : * referred to as `spectral_coefs` in the methods below. For this C++ interface,
107 : * they are saved in a single `DataVector`. To help you index the coefficients
108 : * as expected by this interface, use the class `SpherepackIterator`.
109 : *
110 : * -# nodal: The values at certain collocation points, referred to as
111 : * `collocation_values` in the methods below. This is an array of the expanded
112 : * function \f$g(\theta,\phi)\f$ evaluated at collocation values
113 : * \f$(\theta_i,\phi_j)\f$, where \f$\theta_i\f$ are Gauss-Legendre quadrature
114 : * nodes in the interval \f$(0, \pi)\f$ with \f$i = 0, ..., l_{\max}\f$, and
115 : * \f$\phi_j\f$ is distributed uniformly in \f$(0, 2\pi)\f$ with \f$i = 0, ...,
116 : * 2m_{\max}\f$. The angles of the collocation points can be computed with the
117 : * method `theta_phi_points`.
118 : *
119 : * To convert between the two representations the methods `spec_to_phys` and
120 : * `phys_to_spec` can be used. For internal calculations SPHEREPACK will usually
121 : * convert to spectral coefficients first, so it is in general more efficient to
122 : * use these directly.
123 : *
124 : * Most methods of SPHEREPACK will compute the requested values of e.g.
125 : * `gradient` or `scalar_laplacian` at the collocation points, effectively
126 : * returning an expansion in nodal form as defined above. To evaluate the
127 : * function at arbitrary angles \f$\theta\f$, \f$\phi\f$, these values have to
128 : * be "interpolated" (i.e. the new expansion evaluated) using `interpolate`.
129 : *
130 : * Spherepack stores two types of quantities:
131 : * 1. storage_, which is filled in the constructor and is always const.
132 : * 2. memory_pool_, which is dynamic and thread_local, and is overwritten
133 : * by various member functions that need temporary storage.
134 : */
135 1 : class Spherepack {
136 : public:
137 : /// Type returned by gradient function.
138 1 : using FirstDeriv = tnsr::i<DataVector, 2, Frame::ElementLogical>;
139 : /// Type returned by second derivative function.
140 1 : using SecondDeriv = tnsr::ij<DataVector, 2, Frame::ElementLogical>;
141 :
142 : /// Struct to hold cached information at a set of target interpolation
143 : /// points.
144 : template <typename T>
145 1 : struct InterpolationInfo {
146 0 : InterpolationInfo(size_t l_max, size_t m_max, const gsl::span<double> pmm,
147 : const std::array<T, 2>& target_points);
148 0 : T cos_theta;
149 : // cos(m*phi)
150 0 : DynamicBuffer<T> cos_m_phi;
151 : // sin(m*phi)
152 0 : DynamicBuffer<T> sin_m_phi;
153 : // pbar_factor[m] = Pbar(m,m)*sin(theta)^m
154 0 : DynamicBuffer<T> pbar_factor;
155 :
156 0 : size_t size() const { return num_points_; }
157 0 : size_t m_max() const { return m_max_; }
158 0 : size_t l_max() const { return l_max_; }
159 :
160 : private:
161 0 : size_t l_max_;
162 0 : size_t m_max_;
163 0 : size_t num_points_;
164 : };
165 :
166 : /// Here l_max and m_max are the largest fully-represented l and m in
167 : /// the Ylm expansion.
168 1 : Spherepack(size_t l_max, size_t m_max);
169 :
170 : /// The number of theta collocation points (Gauss-Legendre in
171 : /// \f$\cos\theta\f$) needed to represent all modes up to
172 : /// \f$l_{\max}\f$ exactly.
173 : ///
174 : /// Use to get the correct sizes of vectors of collocation points and
175 : /// spectral coefficients for a given l_max and m_max. Useful for allocating
176 : /// space without having to create a Spherepack.
177 1 : SPECTRE_ALWAYS_INLINE static constexpr size_t n_theta_points(
178 : const size_t l_max) {
179 : return l_max + 1;
180 : }
181 : /// The number of phi collocation points (uniform in \f$\phi\f$) needed
182 : /// to represent all modes up to \f$m_{\max}\f$ exactly.
183 : ///
184 : /// Use to get the correct sizes of vectors of collocation points and
185 : /// spectral coefficients for a given l_max and m_max. Useful for allocating
186 : /// space without having to create a Spherepack.
187 1 : SPECTRE_ALWAYS_INLINE static constexpr size_t n_phi_points(
188 : const size_t m_max) {
189 : return 2 * m_max + 1;
190 : }
191 : /// \brief The number of grid points used in the nodal (physical)
192 : /// representation.
193 : ///
194 : /// Use to get the correct sizes of vectors of collocation points and
195 : /// spectral coefficients for a given l_max and m_max. Useful for allocating
196 : /// space without having to create a Spherepack.
197 1 : SPECTRE_ALWAYS_INLINE static constexpr size_t physical_size(
198 : const size_t l_max, const size_t m_max) {
199 : return n_theta_points(l_max) * n_phi_points(m_max);
200 : }
201 : /// \brief The number of grid points used in the nodal (physical)
202 : /// representation.
203 : ///
204 : /// \note `spectral_size` is the size of the buffer that holds the
205 : /// coefficients; it is not the number of coefficients (which is
206 : /// \f$m_{\max}^2+(l_{\max}-m_{\max})(2m_{\max}+1)\f$).
207 : /// To simplify its internal indexing, SPHEREPACK uses a buffer with
208 : /// more space than necessary. See SpherepackIterator for
209 : /// how to index the coefficients in the buffer.
210 1 : SPECTRE_ALWAYS_INLINE static constexpr size_t spectral_size(
211 : const size_t l_max, const size_t m_max) {
212 : return 2 * (l_max + 1) * (m_max + 1);
213 : }
214 :
215 : /// @{
216 : /// Sizes in physical and spectral space for this instance.
217 1 : size_t l_max() const { return l_max_; }
218 1 : size_t m_max() const { return m_max_; }
219 1 : size_t physical_size() const { return n_theta_ * n_phi_; }
220 1 : size_t spectral_size() const { return spectral_size_; }
221 : /// @}
222 :
223 0 : std::array<size_t, 2> physical_extents() const {
224 : return {{n_theta_, n_phi_}};
225 : }
226 :
227 : /// @{
228 : /// Collocation points theta and phi.
229 : ///
230 : /// The phi points are uniform in phi, with the first point
231 : /// at phi=0.
232 : ///
233 : /// The theta points are Gauss-Legendre in \f$\cos(\theta)\f$,
234 : /// so there are no points at the poles.
235 1 : SPECTRE_ALWAYS_INLINE const std::vector<double>& theta_points() const {
236 : return storage_.theta;
237 : }
238 1 : SPECTRE_ALWAYS_INLINE const std::vector<double>& phi_points() const {
239 : return storage_.phi;
240 : }
241 1 : std::array<DataVector, 2> theta_phi_points() const;
242 : /// @}
243 :
244 : /// @{
245 : /// Spectral transformations.
246 : /// To act on a slice of the input and output arrays, specify strides
247 : /// and offsets.
248 1 : void phys_to_spec(gsl::not_null<double*> spectral_coefs,
249 : gsl::not_null<const double*> collocation_values,
250 : size_t physical_stride = 1, size_t physical_offset = 0,
251 : size_t spectral_stride = 1,
252 : size_t spectral_offset = 0) const {
253 : phys_to_spec_impl(spectral_coefs, collocation_values, physical_stride,
254 : physical_offset, spectral_stride, spectral_offset, false);
255 : }
256 1 : void spec_to_phys(gsl::not_null<double*> collocation_values,
257 : gsl::not_null<const double*> spectral_coefs,
258 : size_t spectral_stride = 1, size_t spectral_offset = 0,
259 : size_t physical_stride = 1,
260 : size_t physical_offset = 0) const {
261 : spec_to_phys_impl(collocation_values, spectral_coefs, spectral_stride,
262 : spectral_offset, physical_stride, physical_offset, false);
263 : };
264 : /// @}
265 :
266 : /// @{
267 : /// Spectral transformations where `collocation_values` and
268 : /// `spectral_coefs` are assumed to point to 3-dimensional arrays
269 : /// (I1 x S2 topology), and the transformations are done for all
270 : /// 'radial' points at once by internally looping over all values of
271 : /// the offset from zero to `stride`-1 (the physical and spectral
272 : /// strides are equal and are called `stride`).
273 1 : void phys_to_spec_all_offsets(gsl::not_null<double*> spectral_coefs,
274 : gsl::not_null<const double*> collocation_values,
275 : size_t stride) const {
276 : phys_to_spec_impl(spectral_coefs, collocation_values, stride, 0, stride, 0,
277 : true);
278 : }
279 1 : void spec_to_phys_all_offsets(gsl::not_null<double*> collocation_values,
280 : gsl::not_null<const double*> spectral_coefs,
281 : size_t stride) const {
282 : spec_to_phys_impl(collocation_values, spectral_coefs, stride, 0, stride, 0,
283 : true);
284 : };
285 : /// @}
286 :
287 : /// @{
288 : /// Simpler, less general interfaces to `phys_to_spec` and `spec_to_phys`.
289 : /// Acts on a slice of the input and returns a unit-stride result.
290 1 : DataVector phys_to_spec(const DataVector& collocation_values,
291 : size_t physical_stride = 1,
292 : size_t physical_offset = 0) const;
293 1 : DataVector spec_to_phys(const DataVector& spectral_coefs,
294 : size_t spectral_stride = 1,
295 : size_t spectral_offset = 0) const;
296 : /// @}
297 :
298 : /// @{
299 : /// Simpler, less general interfaces to `phys_to_spec_all_offsets`
300 : /// and `spec_to_phys_all_offsets`. Result has the same stride as
301 : /// the input.
302 1 : DataVector phys_to_spec_all_offsets(const DataVector& collocation_values,
303 : size_t stride) const;
304 1 : DataVector spec_to_phys_all_offsets(const DataVector& spectral_coefs,
305 : size_t stride) const;
306 : /// @}
307 :
308 : /// Computes Pfaffian derivative (df/dtheta, csc(theta) df/dphi) at
309 : /// the collocation values.
310 : /// To act on a slice of the input and output arrays, specify stride
311 : /// and offset (assumed to be the same for input and output).
312 1 : void gradient(const std::array<double*, 2>& df,
313 : gsl::not_null<const double*> collocation_values,
314 : size_t physical_stride = 1, size_t physical_offset = 0) const;
315 :
316 : /// Same as `gradient`, but takes the spectral coefficients (rather
317 : /// than collocation values) of the function. This is more
318 : /// efficient if one happens to already have the spectral
319 : /// coefficients.
320 : /// To act on a slice of the input and output arrays, specify strides
321 : /// and offsets.
322 1 : void gradient_from_coefs(const std::array<double*, 2>& df,
323 : gsl::not_null<const double*> spectral_coefs,
324 : size_t spectral_stride = 1,
325 : size_t spectral_offset = 0,
326 : size_t physical_stride = 1,
327 : size_t physical_offset = 0) const {
328 : gradient_from_coefs_impl(df, spectral_coefs, spectral_stride,
329 : spectral_offset, physical_stride, physical_offset,
330 : false);
331 : }
332 :
333 : /// @{
334 : /// Same as `gradient` but pointers are assumed to point to
335 : /// 3-dimensional arrays (I1 x S2 topology), and the gradient is
336 : /// done for all 'radial' points at once by internally looping
337 : /// over all values of the offset from zero to `stride`-1.
338 1 : void gradient_all_offsets(const std::array<double*, 2>& df,
339 : gsl::not_null<const double*> collocation_values,
340 : size_t stride = 1) const;
341 :
342 1 : SPECTRE_ALWAYS_INLINE void gradient_from_coefs_all_offsets(
343 : const std::array<double*, 2>& df,
344 : gsl::not_null<const double*> spectral_coefs, size_t stride = 1) const {
345 : gradient_from_coefs_impl(df, spectral_coefs, stride, 0, stride, 0, true);
346 : }
347 : /// @}
348 :
349 : /// @{
350 : /// Simpler, less general interfaces to `gradient`.
351 : /// Acts on a slice of the input and returns a unit-stride result.
352 1 : FirstDeriv gradient(const DataVector& collocation_values,
353 : size_t physical_stride = 1,
354 : size_t physical_offset = 0) const;
355 1 : FirstDeriv gradient_from_coefs(const DataVector& spectral_coefs,
356 : size_t spectral_stride = 1,
357 : size_t spectral_offset = 0) const;
358 : /// @}
359 :
360 : /// @{
361 : /// Simpler, less general interfaces to `gradient_all_offsets`.
362 : /// Result has the same stride as the input.
363 1 : FirstDeriv gradient_all_offsets(const DataVector& collocation_values,
364 : size_t stride = 1) const;
365 1 : FirstDeriv gradient_from_coefs_all_offsets(const DataVector& spectral_coefs,
366 : size_t stride = 1) const;
367 : /// @}
368 :
369 : /// Computes Laplacian in physical space.
370 : /// To act on a slice of the input and output arrays, specify stride
371 : /// and offset (assumed to be the same for input and output).
372 1 : void scalar_laplacian(gsl::not_null<double*> scalar_laplacian,
373 : gsl::not_null<const double*> collocation_values,
374 : size_t physical_stride = 1,
375 : size_t physical_offset = 0) const;
376 :
377 : /// Same as `scalar_laplacian` above, but the input is the spectral
378 : /// coefficients (rather than collocation values) of the function.
379 : /// This is more efficient if one happens to already have the
380 : /// spectral coefficients.
381 : /// To act on a slice of the input and output arrays, specify strides
382 : /// and offsets.
383 1 : void scalar_laplacian_from_coefs(gsl::not_null<double*> scalar_laplacian,
384 : gsl::not_null<const double*> spectral_coefs,
385 : size_t spectral_stride = 1,
386 : size_t spectral_offset = 0,
387 : size_t physical_stride = 1,
388 : size_t physical_offset = 0) const;
389 :
390 : /// @{
391 : /// Simpler, less general interfaces to `scalar_laplacian`.
392 : /// Acts on a slice of the input and returns a unit-stride result.
393 1 : DataVector scalar_laplacian(const DataVector& collocation_values,
394 : size_t physical_stride = 1,
395 : size_t physical_offset = 0) const;
396 1 : DataVector scalar_laplacian_from_coefs(const DataVector& spectral_coefs,
397 : size_t spectral_stride = 1,
398 : size_t spectral_offset = 0) const;
399 : /// @}
400 :
401 : /// Computes Pfaffian first and second derivative in physical space.
402 : /// The first derivative is \f$df(i) = d_i f\f$, and the
403 : /// second derivative is \f$ddf(i,j) = d_i (d_j f)\f$,
404 : /// where \f$d_0 = d/d\theta\f$ and \f$d_1 = csc(\theta) d/d\phi\f$.
405 : /// ddf is not symmetric.
406 : /// To act on a slice of the input and output arrays, specify stride
407 : /// and offset (assumed to be the same for input and output).
408 1 : void second_derivative(const std::array<double*, 2>& df,
409 : gsl::not_null<SecondDeriv*> ddf,
410 : gsl::not_null<const double*> collocation_values,
411 : size_t physical_stride = 1,
412 : size_t physical_offset = 0) const;
413 :
414 : /// Simpler, less general interface to second_derivative
415 1 : std::pair<FirstDeriv, SecondDeriv> first_and_second_derivative(
416 : const DataVector& collocation_values) const;
417 :
418 : /// Computes the integral over the sphere.
419 1 : SPECTRE_ALWAYS_INLINE double definite_integral(
420 : gsl::not_null<const double*> collocation_values,
421 : size_t physical_stride = 1, size_t physical_offset = 0) const {
422 : // clang-tidy: 'do not use pointer arithmetic'
423 : return ddot_(n_theta_ * n_phi_, storage_.quadrature_weights.data(), 1,
424 : collocation_values.get() + physical_offset, // NOLINT
425 : physical_stride);
426 : }
427 :
428 : /// Returns weights \f$w_i\f$ such that \f$sum_i (c_i w_i)\f$
429 : /// is the definite integral, where \f$c_i\f$ are collocation values
430 : /// at point i.
431 1 : SPECTRE_ALWAYS_INLINE const std::vector<double>& integration_weights() const {
432 : return storage_.quadrature_weights;
433 : }
434 :
435 : /// Adds a constant (i.e. \f$f(\theta,\phi)\f$ += \f$c\f$) to the function
436 : /// given by the spectral coefficients, by modifying the coefficients.
437 1 : SPECTRE_ALWAYS_INLINE static void add_constant(
438 : const gsl::not_null<DataVector*> spectral_coefs, const double c) {
439 : // The factor of sqrt(8) is because of the normalization of
440 : // SPHEREPACK's coefficients.
441 : (*spectral_coefs)[0] += sqrt(8.0) * c;
442 : }
443 :
444 : /// Returns the average of \f$f(\theta,\phi)\f$ over \f$(\theta,\phi)\f$.
445 1 : SPECTRE_ALWAYS_INLINE static double average(
446 : const DataVector& spectral_coefs) {
447 : // The factor of sqrt(8) is because of the normalization of
448 : // SPHEREPACK's coefficients. All other coefficients average to zero.
449 : return spectral_coefs[0] / sqrt(8.0);
450 : }
451 :
452 : /// Sets up the `InterpolationInfo` structure for interpolating onto
453 : /// a set of target \f$(\theta,\phi)\f$ points. Does not depend on
454 : /// the function being interpolated.
455 : template <typename T>
456 1 : InterpolationInfo<T> set_up_interpolation_info(
457 : const std::array<T, 2>& target_points) const;
458 :
459 : /// Interpolates from `collocation_values` onto the points that have
460 : /// been passed into the `set_up_interpolation_info` function.
461 : /// To interpolate a different function on the same spectral grid, there
462 : /// is no need to recompute `interpolation_info`.
463 : /// If you specify stride and offset, acts on a slice of the input values.
464 : /// The output has unit stride.
465 : template <typename T>
466 1 : void interpolate(gsl::not_null<T*> result,
467 : gsl::not_null<const double*> collocation_values,
468 : const InterpolationInfo<T>& interpolation_info,
469 : size_t physical_stride = 1,
470 : size_t physical_offset = 0) const;
471 :
472 : /// Same as `interpolate`, but assumes you have spectral coefficients.
473 : /// This is more efficient if you already have the spectral coefficients
474 : /// available.
475 : /// If you specify stride and offset, acts on a slice of the input coefs.
476 : /// The output has unit stride.
477 : template <typename T, typename R>
478 1 : void interpolate_from_coefs(gsl::not_null<T*> result, const R& spectral_coefs,
479 : const InterpolationInfo<T>& interpolation_info,
480 : size_t spectral_stride = 1,
481 : size_t spectral_offset = 0) const;
482 :
483 : /// Simpler interface to `interpolate`. If you need to call this
484 : /// repeatedly on different `spectral_coefs` or `collocation_values`
485 : /// for the same target points, this is inefficient; instead use
486 : /// `set_up_interpolation_info` and the functions that use
487 : /// `InterpolationInfo`.
488 : template <typename T>
489 1 : T interpolate(const DataVector& collocation_values,
490 : const std::array<T, 2>& target_points) const;
491 : template <typename T>
492 0 : T interpolate_from_coefs(const DataVector& spectral_coefs,
493 : const std::array<T, 2>& target_points) const;
494 :
495 : /// Takes spectral coefficients compatible with a resolution given by
496 : /// `l_max_coefs` and `m_max_coefs` and either prolongs them or restricts
497 : /// them to be compatible with a resolution given by `l_max_target` and
498 : /// `m_max_target`. This is done by truncation (restriction) or padding
499 : /// with zeros (prolongation).
500 : ///
501 : /// With `stride > 1`, `spectral_coefs` (and the result) hold `stride`
502 : /// independent sets of coefficients interleaved with the given stride, i.e.
503 : /// coefficient `c` of set `s` is at index `c * stride + s`. This is the
504 : /// layout produced by `phys_to_spec_all_offsets`, so several sets (e.g.
505 : /// radial points) can be prolonged or restricted in a single call.
506 1 : static DataVector prolong_or_restrict(const DataVector& spectral_coefs,
507 : size_t l_max_coefs, size_t m_max_coefs,
508 : size_t l_max_target,
509 : size_t m_max_target, size_t stride = 1);
510 :
511 : /// Same as the overload above, but writes into the pre-allocated `result`
512 : /// buffer (resizing it only if necessary) to avoid allocating. Useful to
513 : /// reuse a single buffer across repeated calls.
514 1 : static void prolong_or_restrict(gsl::not_null<DataVector*> result,
515 : const DataVector& spectral_coefs,
516 : size_t l_max_coefs, size_t m_max_coefs,
517 : size_t l_max_target, size_t m_max_target,
518 : size_t stride = 1);
519 :
520 : /// Takes spectral coefficients compatible with `*this`, and either
521 : /// prolongs them or restricts them to be compatible with `target`.
522 : /// This is done by truncation (restriction) or padding with zeros
523 : /// (prolongation).
524 1 : DataVector prolong_or_restrict(const DataVector& spectral_coefs,
525 : const Spherepack& target) const;
526 :
527 : private:
528 : // Spectral transformations and gradient.
529 : // If `loop_over_offset` is true, then `collocation_values` and
530 : // `spectral_coefs` are assumed to point to 3-dimensional
531 : // arrays (I1 x S2 topology), and the transformations are done for
532 : // all 'radial' points at once by looping over all values of the
533 : // offset from zero to stride-1. If `loop_over_offset` is true,
534 : // `physical_stride` must equal `spectral_stride`.
535 0 : void phys_to_spec_impl(gsl::not_null<double*> spectral_coefs,
536 : gsl::not_null<const double*> collocation_values,
537 : size_t physical_stride = 1, size_t physical_offset = 0,
538 : size_t spectral_stride = 1, size_t spectral_offset = 0,
539 : bool loop_over_offset = false) const;
540 0 : void spec_to_phys_impl(gsl::not_null<double*> collocation_values,
541 : gsl::not_null<const double*> spectral_coefs,
542 : size_t spectral_stride = 1, size_t spectral_offset = 0,
543 : size_t physical_stride = 1, size_t physical_offset = 0,
544 : bool loop_over_offset = false) const;
545 0 : void gradient_from_coefs_impl(const std::array<double*, 2>& df,
546 : gsl::not_null<const double*> spectral_coefs,
547 : size_t spectral_stride = 1,
548 : size_t spectral_offset = 0,
549 : size_t physical_stride = 1,
550 : size_t physical_offset = 0,
551 : bool loop_over_offset = false) const;
552 0 : void calculate_collocation_points();
553 0 : void calculate_interpolation_data();
554 0 : void fill_scalar_work_arrays();
555 0 : void fill_vector_work_arrays();
556 0 : size_t l_max_, m_max_, n_theta_, n_phi_;
557 0 : size_t spectral_size_;
558 : // memory_pool_ will be shared by multiple instances of
559 : // Spherepack on the same thread. Because these instances are on
560 : // the same thread, member functions of two or more of these
561 : // instances cannot be called simultaneously. Note that member
562 : // functions do not make any assumptions about the contents of
563 : // memory_pool_ on entry, so between calls to member functions it is
564 : // safe to resize objects in memory_pool_ or to overwrite them with
565 : // arbitrary data.
566 0 : static thread_local Spherepack_detail::MemoryPool memory_pool_;
567 0 : Spherepack_detail::ConstStorage storage_;
568 : }; // class Spherepack
569 :
570 0 : bool operator==(const Spherepack& lhs, const Spherepack& rhs);
571 0 : bool operator!=(const Spherepack& lhs, const Spherepack& rhs);
572 :
573 : } // namespace ylm
|
