Spectral

Namespaces
namespace	Spectral
	Functionality associated with a particular choice of basis functions and quadrature for spectral operations.
namespace	Spectral::Swsh
	Namespace for spin-weighted spherical harmonic utilities.
namespace	Spectral::filtering
	Matrices for filtering spectral coefficients.

Classes
class	Spectral::Chebyshev
	A collection of helper functions for Chebyshev polynomials. More...
class	Spectral::Fourier
	A collection of helper functions for using a Fourier series. More...
class	ylm::Spherepack
	Defines the C++ interface to SPHEREPACK. More...
class	ylm::SpherepackIterator
	Iterates over spectral coefficients stored in SPHEREPACK format. More...

Typedefs
template<typename DerivativeTagList , ComplexRepresentation Representation = ComplexRepresentation::Interleaved>
using	Spectral::Swsh::AngularDerivatives = detail::AngularDerivativesImpl< DerivativeTagList, typename detail::unique_derived_from_list< DerivativeTagList >::type, Representation >
	A DataBox mutate-compatible computational struct for computing a set of spin-weighted spherical harmonic derivatives by grouping and batch-computing spin-weighted spherical harmonic transforms. More...

Functions
ConvergenceInfo	PowerMonitors::convergence_rate_and_number_of_pile_up_modes (const DataVector &power_monitor, size_t number_of_filtered_modes=0)
	Returns the convergence rate and the number of pile up modes of a power monitor as a ConvergenceInfo. More...
template<int Spin>
auto	Spectral::Swsh::swsh_buffer (const size_t l_max, const size_t number_of_radial_points)
	Produces a SpinWeighted<ComplexModalVector, Spin> of the appropriate size to be used as a modal buffer for Spectral::Swsh::AngularDerivatives or Spectral::Swsh::angular_derivatives. More...
template<typename DerivativeKindList , ComplexRepresentation Representation = ComplexRepresentation::Interleaved, typename... ArgumentTypes>
void	Spectral::Swsh::angular_derivatives (const size_t l_max, const size_t number_of_radial_points, const ArgumentTypes &... arguments)
	Evaluate all of the spin-weighted derivatives in DerivKindList on input SpinWeighted<ComplexDataVector, Spin> collocation data, returning by pointer. More...
template<typename DerivKind , ComplexRepresentation Representation = ComplexRepresentation::Interleaved, int Spin>
SpinWeighted< ComplexDataVector, Tags::derivative_spin_weight< DerivKind >+Spin >	Spectral::Swsh::angular_derivative (size_t l_max, size_t number_of_radial_points, const SpinWeighted< ComplexDataVector, Spin > &to_differentiate)
	Evaluate the spin-weighted derivative DerivKind on the provided SpinWeighted<ComplexDataVector, Spin> collocation data, returning by value.
template<size_t Dim>
void	to_modal_coefficients (gsl::not_null< ComplexModalVector * > modal_coefficients, const ComplexDataVector &nodal_coefficients, const Mesh< Dim > &mesh)
	Compute the modal coefficients from the nodal coefficients. More...
template<size_t Dim>
void	to_modal_coefficients (gsl::not_null< ModalVector * > modal_coefficients, const DataVector &nodal_coefficients, const Mesh< Dim > &mesh)
	Compute the modal coefficients from the nodal coefficients. More...
template<size_t Dim>
ModalVector	to_modal_coefficients (const DataVector &nodal_coefficients, const Mesh< Dim > &mesh)
	Compute the modal coefficients from the nodal coefficients. More...
template<size_t Dim>
ComplexModalVector	to_modal_coefficients (const ComplexDataVector &nodal_coefficients, const Mesh< Dim > &mesh)
	Compute the modal coefficients from the nodal coefficients. More...
template<size_t Dim>
void	to_nodal_coefficients (gsl::not_null< ComplexDataVector * > nodal_coefficients, const ComplexModalVector &modal_coefficients, const Mesh< Dim > &mesh)
	Compute the nodal coefficients from the modal coefficients. More...
template<size_t Dim>
void	to_nodal_coefficients (gsl::not_null< DataVector * > nodal_coefficients, const ModalVector &modal_coefficients, const Mesh< Dim > &mesh)
	Compute the nodal coefficients from the modal coefficients. More...
template<size_t Dim>
DataVector	to_nodal_coefficients (const ModalVector &modal_coefficients, const Mesh< Dim > &mesh)
	Compute the nodal coefficients from the modal coefficients. More...
template<size_t Dim>
ComplexDataVector	to_nodal_coefficients (const ComplexModalVector &modal_coefficients, const Mesh< Dim > &mesh)
	Compute the nodal coefficients from the modal coefficients. More...
template<typename VectorType , size_t Dim>
void	PowerMonitors::power_monitors (gsl::not_null< std::array< DataVector, Dim > * > result, const VectorType &u, const Mesh< Dim > &mesh)
	Returns array of power monitors in each spatial dimension. More...
template<typename VectorType , size_t Dim>
std::array< DataVector, Dim >	PowerMonitors::power_monitors (const VectorType &u, const Mesh< Dim > &mesh)
	Returns array of power monitors in each spatial dimension. More...
double	PowerMonitors::relative_truncation_error (const DataVector &power_monitor, size_t num_modes_to_use)
	Compute the relative truncation error. More...
template<typename VectorType , size_t Dim>
std::array< double, Dim >	PowerMonitors::absolute_truncation_error (const VectorType &tensor_component, const Mesh< Dim > &mesh)
	Returns an estimate of the absolute truncation error in each dimension. More...
void	ylm::real_spherical_harmonic (gsl::not_null< DataVector * > spherical_harmonic, const DataVector &theta, const DataVector &phi, size_t l, int m)
	Evaluates a real spherical harmonic of order l at the requested angles \(\theta\) and \(\phi\). More...
DataVector	ylm::real_spherical_harmonic (const DataVector &theta, const DataVector &phi, size_t l, int m)
	Evaluates a real spherical harmonic of order l at the requested angles \(\theta\) and \(\phi\). More...

Detailed Description

Things related to spectral transformations.

Typedef Documentation

AngularDerivatives

template<typename DerivativeTagList , ComplexRepresentation Representation = ComplexRepresentation::Interleaved>

using Spectral::Swsh::AngularDerivatives = typedef detail::AngularDerivativesImpl< DerivativeTagList, typename detail::unique_derived_from_list<DerivativeTagList>::type, Representation>

A DataBox mutate-compatible computational struct for computing a set of spin-weighted spherical harmonic derivatives by grouping and batch-computing spin-weighted spherical harmonic transforms.

Details

A derivative is evaluated for each tag in DerivativeTagList. All entries in DerivativeTagList must be the tag Spectral::Swsh::Tags::Derivative<Tag, DerivativeKind> prefixing the Tag to be differentiated, and indicating the spin-weighted derivative DerivativeKind to be taken. A DataBox on which this struct is invoked must contain:

• each of the tags in DerivativeTagList (the results of the computation)
• each of the tags Tag prefixed by Spectral::Swsh::Tags::Derivative in DerivativeTagList (the inputs of the computation).
• each of the tags Spectral::Swsh::Tags::SwshTransform<DerivativeTag> for DerivativeTagin DerivativeTagList (the buffers for the derivative applied to the modes)
• each of the tags Spectral::Swsh::Tags::SwshTransform<Tag> for Tag prefixed by any DerivativeTag in DerivativeTagList (the buffers for the transforms of the input data).

This function optimizes the derivative taking process by clustering like spins of tags, forward-transforming each spin cluster together, applying the factor for the derivative to each modal vector, re-clustering according to the new spin weights (the derivatives alter the spin weights), and finally inverse-transforming in clusters.

Function Documentation

absolute_truncation_error()

template<typename VectorType , size_t Dim>

std::array< double, Dim > PowerMonitors::absolute_truncation_error	(	const VectorType &	tensor_component,
		const Mesh< Dim > &	mesh
	)

Returns an estimate of the absolute truncation error in each dimension.

The estimate of the numerical error is given by

\begin{align*} \mathcal{E}\left[P_k\right] = u_\mathrm{max} \times 10^{- \mathcal{T}[P_k]}, \end{align*}

where \( u_\mathrm{max} = \mathrm{max} |u|\) in the corresponding element and \( \mathcal{T}[P_k] \) is the relative error estimate computed from the power monitors \( P_k \).

Warning

This estimate is intended for visualization purposes only.

angular_derivatives()

template<typename DerivativeKindList , ComplexRepresentation Representation = ComplexRepresentation::Interleaved, typename... ArgumentTypes>

void Spectral::Swsh::angular_derivatives	(	const size_t	l_max,
		const size_t	number_of_radial_points,
		const ArgumentTypes &...	arguments
	)

Evaluate all of the spin-weighted derivatives in DerivKindList on input SpinWeighted<ComplexDataVector, Spin> collocation data, returning by pointer.

Details

This function provides two interfaces, one in which the caller provides the intermediate coefficient buffers needed during the computation of the derivatives, and one in which those buffers are temporarily allocated during the derivative function calls.

For the interface in which the caller does not provide buffers, the arguments must take the following structure (enforced by internal function calls):

• size_t l_max : angular resolution for the spherical representation
• size_t number_of_radial_points : radial resolution (number of consecutive blocks to evaluate derivatives, for each input vector )
• for each DerivKind in DerivKindList, a gsl::not_null<SpinWeighted<ComplexDataVector, Spin + Tags::derivative_spin_weight<DerivKind>>> : the output of the derivative evaluation
• for each DerivKind in DerivKindList, a const SpinWeighted<ComplexDataVector, Spin>& (where the Spin for these arguments matches the corresponding vector from the previous set) : the input to the derivative evaluation.

For the interface in which the caller does provide buffers, the arguments must take the following structure (enforced by internal function calls):

• size_t l_max : angular resolution for the spherical representation
• size_t number_of_radial_points : radial resolution (number of consecutive blocks to evaluate derivatives, for each input vector )
• for each DerivKind in DerivKindList, a gsl::not_null<SpinWeighted<ComplexModalVector, Spin + Tags::derivative_spin_weight<DerivKind>>> : the buffer for the spin-weighted spherical harmonic modes of the derivative quantities.
• for each DerivKind in DerivKindList, a const SpinWeighted<ComplexModalVector, Spin> (where the Spin for these arguments matches the corresponding vector from the previous set) : the buffer for the spin-weighted spherical harmonic modes of the input quantities.
• for each DerivKind in DerivKindList, a gsl::not_null<SpinWeighted<ComplexDataVector, Spin + Tags::derivative_spin_weight<DerivKind>>> : the output of the derivative evaluation
• for each DerivKind in DerivKindList, a const SpinWeighted<ComplexDataVector, Spin> (where the Spin for these arguments matches the corresponding vector from the previous set) : the input to the derivative evaluation.

The function swsh_buffer assists in generating the modal buffers of appropriate size.

convergence_rate_and_number_of_pile_up_modes()

ConvergenceInfo PowerMonitors::convergence_rate_and_number_of_pile_up_modes	(	const DataVector &	power_monitor,
		size_t	number_of_filtered_modes = 0
	)

Returns the convergence rate and the number of pile up modes of a power monitor as a ConvergenceInfo.

Details

Computes the convergence rate of a power monitor as a weighted average of slopes measured using different subsets of spectral modes in the power monitor. Equation (53) of [196] gives the convergence rate \(\mathcal{C}\) in terms of a power monitor \(P_k\) as

\begin{equation} \mathcal{C}(P_k) = -\frac{\sum_{k_1=0}^2\sum_{k_2=\tilde{k}_1}^{\tilde{N}-1} \frac{\mathcal{S}(k_1,k_2)}{\epsilon + \mathcal{E}(k_1,k_2)}}{ \sum_{k_1=0}^2\sum_{k_2=\tilde{k}_1}^{\tilde{N}-1} \frac{1}{\epsilon + \mathcal{E}(k_1,k_2)}}. \end{equation}

Here, \(\mathcal{S}(k_1,k_2)\) is the slope of a linear regression fit of \(\log_{10}(P_k)\) with \(k\) satisfying \(k_1\leq k \leq k_2\), \(\mathcal{E}(k_1,k_2)\) is the error of the slope in that fit, \(\epsilon=\max\left(10^{-3}\max\left(\mathcal{E}(k_1,k_2)\right), 10^{-15}\right)\) is a small number to avoid dividing by zero in the event the fit errors vanish, \(\max\left(\mathcal{E}(k_1,k_2)\right)\) is the maximum fit error of each fit whose slope is included in the summation, \(\tilde{k}_1 = \min\left(k_1+4,\tilde{N}-1\right)\), \(\tilde{N} = N-N_f\), \(N\) is the number of modes in the power monitor, and the highest \(N_f\) modes are filtered. Note that the way \(\epsilon\) is defined is so that it matches SpEC's definition, while also ensuring that it is nonzero even if the error in the slope fit is exactly zero.

Also computes the number of pile up modes in a power monitor. Pile up modes are modes where the power is no longer converging at the overall convergence rate. Following Eq. (56) of [196], the number of pile up modes \(\mathcal{P}\) is defined as

\begin{equation} \mathcal{P}(P_k) = \sum_{j=2}^{\tilde{N}-2} \exp\left[-32\left(\frac{\tilde{\mathcal{C}}_j} {\mathcal{C}(P_k)}\right)^2\right], \end{equation}

where \(\mathcal{C}(P_k)\) is the convergence rate of the power monitor \(P_k\), the local convergence rate \(\tilde{\mathcal{C}}_j\) of mode \(j\) is

\begin{equation} \tilde{\mathcal{C}}_j = -\mathcal{S}(j,\min(\tilde{N}-1,j+4)), \end{equation}

\(\mathcal{S}(k_1,k_2)\) is the slope of a linear regression fit of \(\log_{10}(P_k)\) with \(k\) satisfying \(k_1\leq k \leq k_2\), \(\tilde{N} = N-N_f\), \(N\) is the number of modes in the power monitor, and the highest \(N_f\) modes are filtered. The motivation of this definition is the following: if the local convergence rate \(\tilde{\mathcal{C}}_j\) is comparable to the overall convergence rate \(\mathcal{C}(P_k)\), then the \(j^{\rm th}\) term in the summation becomes \(\approx \exp(-32) \approx 10^{-14}\), while if \(\tilde{\mathcal{C}}_j \ll \mathcal{C}(P_k)\), then the \(j^{\rm th}\) term in the summation is \(\approx \exp(0) = 1\). Note that the coefficient value 32 is chosen to agree with SpEC.

Note

The summation goes up to \(\tilde{N}-2\) so that there is it least one larger unfiltered mode for use in computing the slope. The highest mode used when computing the slope is the highest unfiltered mode, \(\tilde{N}-1\). Mode numbers in the power monitor are zero based. These choices are off by one vs. Eqs. (55) and (56) of [196], because those formulas apparently assume one-based indexing.

Parameters

power_monitor	The power monitor.
number_of_filtered_modes	How many of the highest modes of the power monitor are filtered (default 0).

power_monitors() [1/2]

template<typename VectorType , size_t Dim>

std::array< DataVector, Dim > PowerMonitors::power_monitors	(	const VectorType &	u,
		const Mesh< Dim > &	mesh
	)

Returns array of power monitors in each spatial dimension.

Computed following Sec. 5.1 of Ref. [196]. For example, in the x dimension (indexed by \( k_0 \)), we compute

\begin{align*} P_{k_0}[\psi] = \sqrt{ \frac{1}{N_1 N_2} \sum_{k_1,k_2} \left| C_{k_0,k_1,k_2} \right|^2} , \end{align*}

where \( C_{k_0,k_1,k_2}\) are the modal coefficients of variable \( \psi \).

power_monitors() [2/2]

template<typename VectorType , size_t Dim>

void PowerMonitors::power_monitors	(	gsl::not_null< std::array< DataVector, Dim > * >	result,
		const VectorType &	u,
		const Mesh< Dim > &	mesh
	)

Returns array of power monitors in each spatial dimension.

Computed following Sec. 5.1 of Ref. [196]. For example, in the x dimension (indexed by \( k_0 \)), we compute

\begin{align*} P_{k_0}[\psi] = \sqrt{ \frac{1}{N_1 N_2} \sum_{k_1,k_2} \left| C_{k_0,k_1,k_2} \right|^2} , \end{align*}

where \( C_{k_0,k_1,k_2}\) are the modal coefficients of variable \( \psi \).

real_spherical_harmonic() [1/2]

DataVector ylm::real_spherical_harmonic	(	const DataVector &	theta,
		const DataVector &	phi,
		size_t	l,
		int	m
	)

Evaluates a real spherical harmonic of order l at the requested angles \(\theta\) and \(\phi\).

The real spherical harmonics are defined as:

\begin{equation} Y_{lm}(\theta, \phi) = \begin{cases} \sqrt{2} (-1)^m \mathcal{I}[Y_l^{|m|}] & \text{if } m < 0 \\ Y_l^{0} & \text{if } m = 0 \\ \sqrt{2} (-1)^m \mathcal{R}[Y_l^{m}] & \text{if } m > 0 \end{cases} \end{equation}

where \(Y_l^m\) are the complex spherical harmonics and \(\mathcal{R}\) denotes the real part and \(\mathcal{I}\) denotes the imaginary part.

Note

This implementation uses the boost implementation of spherical harmonics. The calculation is not vectorized and may be slow. Stability has not been tested for large l.

real_spherical_harmonic() [2/2]

void ylm::real_spherical_harmonic	(	gsl::not_null< DataVector * >	spherical_harmonic,
		const DataVector &	theta,
		const DataVector &	phi,
		size_t	l,
		int	m
	)

Evaluates a real spherical harmonic of order l at the requested angles \(\theta\) and \(\phi\).

The real spherical harmonics are defined as:

\begin{equation} Y_{lm}(\theta, \phi) = \begin{cases} \sqrt{2} (-1)^m \mathcal{I}[Y_l^{|m|}] & \text{if } m < 0 \\ Y_l^{0} & \text{if } m = 0 \\ \sqrt{2} (-1)^m \mathcal{R}[Y_l^{m}] & \text{if } m > 0 \end{cases} \end{equation}

where \(Y_l^m\) are the complex spherical harmonics and \(\mathcal{R}\) denotes the real part and \(\mathcal{I}\) denotes the imaginary part.

Note

This implementation uses the boost implementation of spherical harmonics. The calculation is not vectorized and may be slow. Stability has not been tested for large l.

relative_truncation_error()

double PowerMonitors::relative_truncation_error	(	const DataVector &	power_monitor,
		size_t	num_modes_to_use
	)

Compute the relative truncation error.

The negative logarithm of this quantity is defined by Eqs. (57) and (58) of Ref. [196], i.e.,

\begin{align*} \mathcal{T}\left[P_k\right] = \log_{10} \max \left(P_0, P_1\right) - \dfrac{\sum_{j=0}^{j_{\text{max}, k}} \log_{10} \left(P_j\right) w_j} {\sum_{j=0}^{j_{\text{max}, k}} w_j} , \end{align*}

with weights

\begin{align*} w_j = \exp\left[ - \left(j - j_{\text{max}, k} + \dfrac{1}{2}\right)^2 \right] . \end{align*}

where \( j_{\text{max}, k} = N_k - 1 \) and \( N_k \) is the number of modes or gridpoints in dimension k. Here the second term is a weighted average with larger weights toward the highest modes.

Note

Modes below a cutoff of \(100 \epsilon \mathrm{max}_k(P_k)\) are ignored in the weighted average, where \(\epsilon\) is the machine epsilon. This ensures that we don't underestimate the truncation error if some modes are zero (e.g. by symmetry). Furthermore, if the last two or more modes are zero, we assume that the function is represented exactly and return a relative truncation error of zero.

Details

The number of modes (num_modes_to_use) argument needs to be less or equal than the total number of power monitors (power_monitor.size()). In contrast with Ref. [196], here we index the modes starting from zero.

swsh_buffer()

template<int Spin>

auto Spectral::Swsh::swsh_buffer	(	const size_t	l_max,
		const size_t	number_of_radial_points
	)

Produces a SpinWeighted<ComplexModalVector, Spin> of the appropriate size to be used as a modal buffer for Spectral::Swsh::AngularDerivatives or Spectral::Swsh::angular_derivatives.

Details

The Spectral::Swsh::angular_derivatives and Spectral::Swsh::AngularDerivatives interfaces require that calling code provides a buffer for the intermediate transform results, to ensure that callers are aware of the allocations and can suitably reuse buffers if possible. This utility eases the creation of those buffers.

to_modal_coefficients() [1/4]

template<size_t Dim>

ComplexModalVector to_modal_coefficients	(	const ComplexDataVector &	nodal_coefficients,
		const Mesh< Dim > &	mesh
	)

Compute the modal coefficients from the nodal coefficients.

See also

Spectral::nodal_to_modal_matrix

to_modal_coefficients() [2/4]

template<size_t Dim>

ModalVector to_modal_coefficients	(	const DataVector &	nodal_coefficients,
		const Mesh< Dim > &	mesh
	)

Compute the modal coefficients from the nodal coefficients.

See also

Spectral::nodal_to_modal_matrix

to_modal_coefficients() [3/4]

template<size_t Dim>

void to_modal_coefficients	(	gsl::not_null< ComplexModalVector * >	modal_coefficients,
		const ComplexDataVector &	nodal_coefficients,
		const Mesh< Dim > &	mesh
	)

Compute the modal coefficients from the nodal coefficients.

See also

Spectral::nodal_to_modal_matrix

to_modal_coefficients() [4/4]

template<size_t Dim>

void to_modal_coefficients	(	gsl::not_null< ModalVector * >	modal_coefficients,
		const DataVector &	nodal_coefficients,
		const Mesh< Dim > &	mesh
	)

Compute the modal coefficients from the nodal coefficients.

See also

Spectral::nodal_to_modal_matrix

to_nodal_coefficients() [1/4]

template<size_t Dim>

ComplexDataVector to_nodal_coefficients	(	const ComplexModalVector &	modal_coefficients,
		const Mesh< Dim > &	mesh
	)

Compute the nodal coefficients from the modal coefficients.

See also

Spectral::modal_to_nodal_matrix

to_nodal_coefficients() [2/4]

template<size_t Dim>

DataVector to_nodal_coefficients	(	const ModalVector &	modal_coefficients,
		const Mesh< Dim > &	mesh
	)

Compute the nodal coefficients from the modal coefficients.

See also

Spectral::modal_to_nodal_matrix

to_nodal_coefficients() [3/4]

template<size_t Dim>

void to_nodal_coefficients	(	gsl::not_null< ComplexDataVector * >	nodal_coefficients,
		const ComplexModalVector &	modal_coefficients,
		const Mesh< Dim > &	mesh
	)

Compute the nodal coefficients from the modal coefficients.

See also

Spectral::modal_to_nodal_matrix

to_nodal_coefficients() [4/4]

template<size_t Dim>

void to_nodal_coefficients	(	gsl::not_null< DataVector * >	nodal_coefficients,
		const ModalVector &	modal_coefficients,
		const Mesh< Dim > &	mesh
	)

Compute the nodal coefficients from the modal coefficients.

See also

Spectral::modal_to_nodal_matrix