|
SpECTRE
v2024.05.11
|
Defines the C++ interface to SPHEREPACK. More...
#include <Spherepack.hpp>
Classes | |
| struct | InterpolationInfo |
| Struct to hold cached information at a set of target interpolation points. More... | |
Public Types | |
| using | FirstDeriv = tnsr::i< DataVector, 2, Frame::ElementLogical > |
| Type returned by gradient function. | |
| using | SecondDeriv = tnsr::ij< DataVector, 2, Frame::ElementLogical > |
| Type returned by second derivative function. | |
Public Member Functions | |
| Spherepack (size_t l_max, size_t m_max) | |
| Here l_max and m_max are the largest fully-represented l and m in the Ylm expansion. | |
| std::array< size_t, 2 > | physical_extents () const |
| 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 |
| 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). | |
| 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 |
Same as gradient, but takes the spectral coefficients (rather than collocation values) of the function. This is more efficient if one happens to already have the spectral coefficients. To act on a slice of the input and output arrays, specify strides and offsets. | |
| 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 |
| 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). | |
| 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 |
Same as scalar_laplacian above, but the input is the spectral coefficients (rather than collocation values) of the function. This is more efficient if one happens to already have the spectral coefficients. To act on a slice of the input and output arrays, specify strides and offsets. | |
| 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 |
| Computes Pfaffian first and second derivative in physical space. The first derivative is \(df(i) = d_i f\), and the second derivative is \(ddf(i,j) = d_i (d_j f)\), where \(d_0 = d/d\theta\) and \(d_1 = csc(\theta) d/d\phi\). ddf is not symmetric. To act on a slice of the input and output arrays, specify stride and offset (assumed to be the same for input and output). | |
| std::pair< FirstDeriv, SecondDeriv > | first_and_second_derivative (const DataVector &collocation_values) const |
| Simpler, less general interface to second_derivative. | |
| double | definite_integral (gsl::not_null< const double * > collocation_values, size_t physical_stride=1, size_t physical_offset=0) const |
| Computes the integral over the sphere. | |
| const std::vector< double > & | integration_weights () const |
| Returns weights \(w_i\) such that \(sum_i (c_i w_i)\) is the definite integral, where \(c_i\) are collocation values at point i. | |
| template<typename T > | |
| InterpolationInfo< T > | set_up_interpolation_info (const std::array< T, 2 > &target_points) const |
Sets up the InterpolationInfo structure for interpolating onto a set of target \((\theta,\phi)\) points. Does not depend on the function being interpolated. | |
| template<typename T > | |
| void | interpolate (gsl::not_null< T * > result, gsl::not_null< const double * > collocation_values, const InterpolationInfo< T > &interpolation_info, size_t physical_stride=1, size_t physical_offset=0) const |
Interpolates from collocation_values onto the points that have been passed into the set_up_interpolation_info function. To interpolate a different function on the same spectral grid, there is no need to recompute interpolation_info. If you specify stride and offset, acts on a slice of the input values. The output has unit stride. | |
| template<typename T , typename R > | |
| void | interpolate_from_coefs (gsl::not_null< T * > result, const R &spectral_coefs, const InterpolationInfo< T > &interpolation_info, size_t spectral_stride=1, size_t spectral_offset=0) const |
Same as interpolate, but assumes you have spectral coefficients. This is more efficient if you already have the spectral coefficients available. If you specify stride and offset, acts on a slice of the input coefs. The output has unit stride. | |
| template<typename T > | |
| T | interpolate (const DataVector &collocation_values, const std::array< T, 2 > &target_points) const |
Simpler interface to interpolate. If you need to call this repeatedly on different spectral_coefs or collocation_values for the same target points, this is inefficient; instead use set_up_interpolation_info and the functions that use InterpolationInfo. | |
| template<typename T > | |
| T | interpolate_from_coefs (const DataVector &spectral_coefs, const std::array< T, 2 > &target_points) const |
| DataVector | prolong_or_restrict (const DataVector &spectral_coefs, const Spherepack &target) const |
Takes spectral coefficients compatible with *this, and either prolongs them or restricts them to be compatible with target. This is done by truncation (restriction) or padding with zeros (prolongation). | |
| size_t | l_max () const |
| Sizes in physical and spectral space for this instance. | |
| size_t | m_max () const |
| Sizes in physical and spectral space for this instance. | |
| size_t | physical_size () const |
| Sizes in physical and spectral space for this instance. | |
| size_t | spectral_size () const |
| Sizes in physical and spectral space for this instance. | |
| const std::vector< double > & | theta_points () const |
| Collocation points theta and phi. More... | |
| const std::vector< double > & | phi_points () const |
| Collocation points theta and phi. More... | |
| std::array< DataVector, 2 > | theta_phi_points () const |
| Collocation points theta and phi. More... | |
| 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 |
| Spectral transformations. To act on a slice of the input and output arrays, specify strides and offsets. | |
| 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 |
| Spectral transformations. To act on a slice of the input and output arrays, specify strides and offsets. | |
| void | phys_to_spec_all_offsets (gsl::not_null< double * > spectral_coefs, gsl::not_null< const double * > collocation_values, size_t stride) const |
Spectral transformations where collocation_values and spectral_coefs are assumed to point to 3-dimensional arrays (I1 x S2 topology), and the transformations are done for all 'radial' points at once by internally looping over all values of the offset from zero to stride-1 (the physical and spectral strides are equal and are called stride). | |
| void | spec_to_phys_all_offsets (gsl::not_null< double * > collocation_values, gsl::not_null< const double * > spectral_coefs, size_t stride) const |
Spectral transformations where collocation_values and spectral_coefs are assumed to point to 3-dimensional arrays (I1 x S2 topology), and the transformations are done for all 'radial' points at once by internally looping over all values of the offset from zero to stride-1 (the physical and spectral strides are equal and are called stride). | |
| DataVector | phys_to_spec (const DataVector &collocation_values, size_t physical_stride=1, size_t physical_offset=0) const |
Simpler, less general interfaces to phys_to_spec and spec_to_phys. Acts on a slice of the input and returns a unit-stride result. | |
| DataVector | spec_to_phys (const DataVector &spectral_coefs, size_t spectral_stride=1, size_t spectral_offset=0) const |
Simpler, less general interfaces to phys_to_spec and spec_to_phys. Acts on a slice of the input and returns a unit-stride result. | |
| DataVector | phys_to_spec_all_offsets (const DataVector &collocation_values, size_t stride) const |
Simpler, less general interfaces to phys_to_spec_all_offsets and spec_to_phys_all_offsets. Result has the same stride as the input. | |
| DataVector | spec_to_phys_all_offsets (const DataVector &spectral_coefs, size_t stride) const |
Simpler, less general interfaces to phys_to_spec_all_offsets and spec_to_phys_all_offsets. Result has the same stride as the input. | |
| void | gradient_all_offsets (const std::array< double *, 2 > &df, gsl::not_null< const double * > collocation_values, size_t stride=1) const |
Same as gradient but pointers are assumed to point to 3-dimensional arrays (I1 x S2 topology), and the gradient is done for all 'radial' points at once by internally looping over all values of the offset from zero to stride-1. | |
| void | gradient_from_coefs_all_offsets (const std::array< double *, 2 > &df, gsl::not_null< const double * > spectral_coefs, size_t stride=1) const |
Same as gradient but pointers are assumed to point to 3-dimensional arrays (I1 x S2 topology), and the gradient is done for all 'radial' points at once by internally looping over all values of the offset from zero to stride-1. | |
| FirstDeriv | gradient (const DataVector &collocation_values, size_t physical_stride=1, size_t physical_offset=0) const |
Simpler, less general interfaces to gradient. Acts on a slice of the input and returns a unit-stride result. | |
| FirstDeriv | gradient_from_coefs (const DataVector &spectral_coefs, size_t spectral_stride=1, size_t spectral_offset=0) const |
Simpler, less general interfaces to gradient. Acts on a slice of the input and returns a unit-stride result. | |
| FirstDeriv | gradient_all_offsets (const DataVector &collocation_values, size_t stride=1) const |
Simpler, less general interfaces to gradient_all_offsets. Result has the same stride as the input. | |
| FirstDeriv | gradient_from_coefs_all_offsets (const DataVector &spectral_coefs, size_t stride=1) const |
Simpler, less general interfaces to gradient_all_offsets. Result has the same stride as the input. | |
| DataVector | scalar_laplacian (const DataVector &collocation_values, size_t physical_stride=1, size_t physical_offset=0) const |
Simpler, less general interfaces to scalar_laplacian. Acts on a slice of the input and returns a unit-stride result. | |
| DataVector | scalar_laplacian_from_coefs (const DataVector &spectral_coefs, size_t spectral_stride=1, size_t spectral_offset=0) const |
Simpler, less general interfaces to scalar_laplacian. Acts on a slice of the input and returns a unit-stride result. | |
Static Public Member Functions | |
| static void | add_constant (const gsl::not_null< DataVector * > spectral_coefs, const double c) |
| Adds a constant (i.e. \(f(\theta,\phi)\) += \(c\)) to the function given by the spectral coefficients, by modifying the coefficients. | |
| static double | average (const DataVector &spectral_coefs) |
| Returns the average of \(f(\theta,\phi)\) over \((\theta,\phi)\). | |
| static constexpr size_t | physical_size (const size_t l_max, const size_t m_max) |
| Static functions to return the correct sizes of vectors of collocation points and spectral coefficients for a given l_max and m_max. Useful for allocating space without having to create a Spherepack. | |
| static constexpr size_t | spectral_size (const size_t l_max, const size_t m_max) |
Defines the C++ interface to SPHEREPACK.
The class Spherepack defines the C++ interface to the fortran library SPHEREPACK used for computations on the surface of a sphere.
Given a real-valued, scalar function \(g(\theta, \phi)\), SPHEREPACK expands it as:
\begin{align} g(\theta, \phi) &=\frac{1}{2}\sum_{l=0}^{l_{\max}}\bar P_l^0(\cos\theta) a_{l0} +\sum_{l=1}^{l_{\max}}\sum_{m=1}^{\min(l, m_{\max})}\bar P_l^m(\cos\theta)\{ a_{lm}\cos m\phi -b_{lm}\sin m\phi\}\label{eq:spherepack_expansion} \end{align}
where \(a_{lm}\) and \(b_{lm}\) are real-valued spectral coefficient arrays used by SPHEREPACK, \(P_l^m(x)\) are defined as
\begin{align} \bar P_l^m(x)&=\sqrt{\frac{(2l+1)(l-m)!}{2(l+m)!}}\;P_{lm}(x) \end{align}
and \(P_{nm}(x)\) are the associated Legendre polynomials as defined, for example, in Jackson's "Classical Electrodynamics".
The standard expansion of \(g(\theta, \phi)\) in terms of scalar spherical harmonics is
\begin{align} g(\theta, \phi) &= \sum_{l=0}^{l_{\max}}\sum_{m=-\min(l, m_{\max})}^{\min(l, m_{\max})} A_{lm} Y_{lm}(\theta,\phi), \end{align}
where \(Y_{lm}(\theta,\phi)\) are the usual complex-valued scalar spherical harmonics (as defined, for example, in Jackson's "Classical Electrodynamics") and \(A_{lm}\) are complex coefficients.
The relationship between the complex coefficients \(A_{lm}\) and SPHEREPACK's real-valued \(a_{lm}\) and \(b_{lm}\) is
\begin{align} a_{l0} & = \sqrt{\frac{2}{\pi}}A_{l0}&\qquad l\geq 0,\\ a_{lm} & = (-1)^m\sqrt{\frac{2}{\pi}} \mathrm{Re}(A_{lm}) &\qquad l\geq 1, m\geq 1, \\ b_{lm} & = (-1)^m\sqrt{\frac{2}{\pi}} \mathrm{Im}(A_{lm}) &\qquad l\geq 1, m\geq 1. \end{align}
Sometimes it is useful to expand a real-valued function in the form
\begin{align} g(\theta, \phi) &= \sum_{l=0}^\infty\sum_{m=0}^l \left[ c_{lm}\mathrm{Re}(Y_{lm}(\theta, \phi))+ d_{nm}\mathrm{Im}(Y_{lm}(\theta, \phi)) \right]. \end{align}
The coefficients here are therefore
\begin{align} c_{l0} &= A_{l0},\\ c_{lm} &= 2\mathrm{Re}(A_{lm}) \qquad m\geq 1,\\ d_{lm} &=-2\mathrm{Im}(A_{lm}). \end{align}
Internally, SPHEREPACK can represent its expansion in two ways which we will refer to as modal and nodal representations:
spectral_coefs in the methods below. For this C++ interface, they are saved in a single DataVector. To help you index the coefficients as expected by this interface, use the class SpherepackIterator.collocation_values in the methods below. This is an array of the expanded function \(g(\theta,\phi)\) evaluated at collocation values \((\theta_i,\phi_j)\), where \(\theta_i\) are Gauss-Legendre quadrature nodes in the interval \((0, \pi)\) with \(i = 0, ..., l_{\max}\), and \(\phi_j\) is distributed uniformly in \((0, 2\pi)\) with \(i = 0, ..., 2m_{\max}\). The angles of the collocation points can be computed with the method theta_phi_points.To convert between the two representations the methods spec_to_phys and phys_to_spec can be used. For internal calculations SPHEREPACK will usually convert to spectral coefficients first, so it is in general more efficient to use these directly.
Most methods of SPHEREPACK will compute the requested values of e.g. gradient or scalar_laplacian at the collocation points, effectively returning an expansion in nodal form as defined above. To evaluate the function at arbitrary angles \(\theta\), \(\phi\), these values have to be "interpolated" (i.e. the new expansion evaluated) using interpolate.
Spherepack stores two types of quantities:
|
inline |
Collocation points theta and phi.
The phi points are uniform in phi, with the first point at phi=0.
The theta points are Gauss-Legendre in \(\cos(\theta)\), so there are no points at the poles.
|
inlinestaticconstexpr |
spectral_size is the size of the buffer that holds the coefficients; it is not the number of coefficients (which is \(m_{\rm max}^2+(l_{\rm max}-m_{\rm max})(2m_{\rm max}+1)\)). To simplify its internal indexing, SPHEREPACK uses a buffer with more space than necessary. See SpherepackIterator for how to index the coefficients in the buffer. | std::array< DataVector, 2 > ylm::Spherepack::theta_phi_points | ( | ) | const |
Collocation points theta and phi.
The phi points are uniform in phi, with the first point at phi=0.
The theta points are Gauss-Legendre in \(\cos(\theta)\), so there are no points at the poles.
|
inline |
Collocation points theta and phi.
The phi points are uniform in phi, with the first point at phi=0.
The theta points are Gauss-Legendre in \(\cos(\theta)\), so there are no points at the poles.