SpECTRE  v2024.06.05
Spectral Namespace Reference

Functionality associated with a particular choice of basis functions and quadrature for spectral operations. More...

Namespaces

namespace  filtering
 Matrices for filtering spectral coefficients.
 
namespace  Swsh
 Namespace for spin-weighted spherical harmonic utilities.
 

Typedefs

using MortarSize = ChildSize
 The portion of an element covered by a mortar.
 

Enumerations

enum class  Basis : uint8_t {
  Uninitialized = 0 << basis_shift , Chebyshev = 1 << basis_shift , Legendre = 2 << basis_shift , FiniteDifference = 3 << basis_shift ,
  SphericalHarmonic = 4 << basis_shift
}
 Either the basis functions used by a spectral or discontinuous Galerkin (DG) method, or the value FiniteDifference when a finite difference method is used. More...
 
enum class  ChildSize : uint8_t { Uninitialized = 0 , Full = 1 , UpperHalf = 2 , LowerHalf = 3 }
 The portion of a mesh covered by a child mesh.
 
enum class  Quadrature : uint8_t {
  Uninitialized = 0 , Gauss , GaussLobatto , CellCentered ,
  FaceCentered , Equiangular
}
 Either the choice of quadrature method to compute integration weights for a spectral or discontinuous Galerkin (DG) method, or the locations of grid points when a finite difference method is used. More...
 

Functions

std::array< Basis, 5 > all_bases ()
 All possible values of Basis.
 
Basis to_basis (const std::string &basis)
 Convert a string to a Basis enum.
 
std::ostreamoperator<< (std::ostream &os, const Basis &basis)
 Output operator for a Basis.
 
template<typename CoefficientDataType , typename DataType >
DataType evaluate_clenshaw (const std::vector< CoefficientDataType > &coefficients, const DataType &alpha, const DataType &beta, const DataType &f0_of_x, const DataType &f1_of_x)
 Clenshaw's algorithm for evaluating linear combinations of basis functions obeying a given three term recurrence relation. See numerical recipes 3rd edition 5.4.2 for details. alpha and beta define the recursion \( f_n(x) = \alpha f_{n-1}(x) + \beta f_{n-2}(x) \) f0_of_x and f1_of_x are \( f_0(x)\) and \(f_1(x)\) respectively.
 
std::ostreamoperator<< (std::ostream &os, ChildSize mortar_size)
 
template<size_t Dim>
bool needs_projection (const Mesh< Dim > &mesh1, const Mesh< Dim > &mesh2, const std::array< ChildSize, Dim > &child_sizes)
 Determine whether data needs to be projected between a child mesh and its parent mesh. If no projection is necessary the data may be used as-is. Projection is necessary if the child is either p-refined or h-refined relative to its parent, or both. This operation is symmetric, i.e. it is irrelevant in which order the child and the parent mesh are passed in.
 
const Matrixprojection_matrix_child_to_parent (const Mesh< 1 > &child_mesh, const Mesh< 1 > &parent_mesh, ChildSize size, bool operand_is_massive=false)
 The projection matrix from a child mesh to its parent. More...
 
template<size_t Dim>
std::array< std::reference_wrapper< const Matrix >, Dim > projection_matrix_child_to_parent (const Mesh< Dim > &child_mesh, const Mesh< Dim > &parent_mesh, const std::array< ChildSize, Dim > &child_sizes, bool operand_is_massive=false)
 The projection matrix from a child mesh to its parent, in Dim dimensions.
 
const Matrixprojection_matrix_parent_to_child (const Mesh< 1 > &parent_mesh, const Mesh< 1 > &child_mesh, ChildSize size)
 The projection matrix from a parent mesh to one of its children. More...
 
template<size_t Dim>
std::array< std::reference_wrapper< const Matrix >, Dim > projection_matrix_parent_to_child (const Mesh< Dim > &parent_mesh, const Mesh< Dim > &child_mesh, const std::array< ChildSize, Dim > &child_sizes)
 The projection matrix from a parent mesh to one of its children, in Dim dimensions.
 
template<size_t Dim>
std::array< std::reference_wrapper< const Matrix >, Dim > p_projection_matrices (const Mesh< Dim > &source_mesh, const Mesh< Dim > &target_mesh)
 The projection matrices from a source mesh to a target mesh where the meshes cover the same physical volume.
 
std::array< Quadrature, 6 > all_quadratures ()
 All possible values of Quadrature.
 
Quadrature to_quadrature (const std::string &quadrature)
 Convert a string to a Quadrature enum.
 
std::ostreamoperator<< (std::ostream &os, const Quadrature &quadrature)
 Output operator for a Quadrature.
 
template<Basis BasisType, typename T >
compute_basis_function_value (size_t k, const T &x)
 Compute the function values of the basis function \(\Phi_k(x)\) (zero-indexed).
 
template<Basis >
DataVector compute_inverse_weight_function_values (const DataVector &)
 Compute the inverse of the weight function \(w(x)\) w.r.t. which the basis functions are orthogonal. See the description of quadrature_weights(size_t) for details. This is arbitrarily set to 1 for FiniteDifference basis, to integrate using the midpoint method (see quadrature_weights (size_t) for details).
 
template<Basis BasisType>
double compute_basis_function_normalization_square (size_t k)
 Compute the normalization square of the basis function \(\Phi_k\) (zero-indexed), i.e. the weighted definite integral over its square.
 
template<Basis BasisType, Quadrature QuadratureType>
std::pair< DataVector, DataVectorcompute_collocation_points_and_weights (size_t num_points)
 Compute the collocation points and weights associated to the basis and quadrature. More...
 
template<Basis BasisType, Quadrature QuadratureType>
const DataVectorcollocation_points (size_t num_points)
 Collocation points. More...
 
const DataVectorcollocation_points (const Mesh< 1 > &mesh)
 Collocation points for a one-dimensional mesh. More...
 
template<Basis BasisType, Quadrature QuadratureType>
const DataVectorquadrature_weights (size_t num_points)
 Weights to compute definite integrals. More...
 
const DataVectorquadrature_weights (const Mesh< 1 > &mesh)
 Quadrature weights for a one-dimensional mesh. More...
 
template<Basis BasisType, Quadrature QuadratureType>
const Matrixweak_flux_differentiation_matrix (size_t num_points)
 Matrix used to compute the divergence of the flux in weak form. More...
 
const Matrixweak_flux_differentiation_matrix (const Mesh< 1 > &mesh)
 Matrix used to compute the divergence of the flux in weak form. More...
 
template<Basis BasisType, Quadrature QuadratureType>
const Matrixintegration_matrix (size_t num_points)
 Matrix used to perform an indefinite integral of a function over the logical grid. The left boundary condition is such that the integral is 0 at \(\xi=-1\) More...
 
const Matrixintegration_matrix (const Mesh< 1 > &mesh)
 Indefinite integration matrix for a one-dimensional mesh. More...
 
template<Basis BasisType, Quadrature QuadratureType, typename T >
Matrix interpolation_matrix (size_t num_points, const T &target_points)
 Matrix used to interpolate to the target_points. More...
 
template<typename T >
Matrix interpolation_matrix (const Mesh< 1 > &mesh, const T &target_points)
 Interpolation matrix to the target_points for a one-dimensional mesh. More...
 
template<Basis BasisType, Quadrature QuadratureType>
const Matrixmodal_to_nodal_matrix (size_t num_points)
 Matrix used to transform from the spectral coefficients (modes) of a function to its nodal coefficients. Also referred to as the Vandermonde matrix. More...
 
const Matrixmodal_to_nodal_matrix (const Mesh< 1 > &mesh)
 Transformation matrix from modal to nodal coefficients for a one-dimensional mesh. More...
 
template<Basis BasisType, Quadrature QuadratureType>
const Matrixnodal_to_modal_matrix (size_t num_points)
 Matrix used to transform from the nodal coefficients of a function to its spectral coefficients (modes). Also referred to as the inverse Vandermonde matrix. More...
 
const Matrixnodal_to_modal_matrix (const Mesh< 1 > &mesh)
 Transformation matrix from nodal to modal coefficients for a one-dimensional mesh. More...
 
template<Basis BasisType, Quadrature QuadratureType>
const Matrixlinear_filter_matrix (size_t num_points)
 Matrix used to linearize a function. More...
 
const Matrixlinear_filter_matrix (const Mesh< 1 > &mesh)
 Linear filter matrix for a one-dimensional mesh. More...
 
template<Basis BasisType, Quadrature QuadratureType>
const Matrixdifferentiation_matrix (size_t num_points)
 Matrix used to compute the derivative of a function. More...
 
template<Basis BasisType, Quadrature QuadratureType>
const Matrixdifferentiation_matrix_transpose (size_t num_points)
 Matrix used to compute the derivative of a function. More...
 
const Matrixdifferentiation_matrix (const Mesh< 1 > &mesh)
 Differentiation matrix for a one-dimensional mesh. More...
 
const Matrixdifferentiation_matrix_transpose (const Mesh< 1 > &mesh)
 Differentiation matrix for a one-dimensional mesh. More...
 
const std::pair< Matrix, Matrix > & boundary_interpolation_matrices (const Mesh< 1 > &mesh)
 Matrices that interpolate to the lower and upper boundaries of the element. More...
 
template<Basis BasisType, Quadrature QuadratureType>
const std::pair< Matrix, Matrix > & boundary_interpolation_matrices (size_t num_points)
 Matrices that interpolate to the lower and upper boundaries of the element. More...
 
const std::pair< DataVector, DataVector > & boundary_interpolation_term (const Mesh< 1 > &mesh)
 Interpolates values from the boundary into the volume, which is needed when applying time derivative or Bjorhus-type boundary conditions in a discontinuous Galerkin scheme using Gauss points. More...
 
template<Basis BasisType, Quadrature QuadratureType>
const std::pair< DataVector, DataVector > & boundary_interpolation_term (size_t num_points)
 Interpolates values from the boundary into the volume, which is needed when applying time derivative or Bjorhus-type boundary conditions in a discontinuous Galerkin scheme using Gauss points. More...
 
const std::pair< DataVector, DataVector > & boundary_lifting_term (const Mesh< 1 > &mesh)
 Terms used during the lifting portion of a discontinuous Galerkin scheme when using Gauss points. More...
 
template<Basis BasisType, Quadrature QuadratureType>
const std::pair< DataVector, DataVector > & boundary_lifting_term (size_t num_points)
 Terms used during the lifting portion of a discontinuous Galerkin scheme when using Gauss points. More...
 

Variables

constexpr uint8_t basis_shift = 4
 The amount we left-shift the basis integers by. More...
 
template<Basis basis, Quadrature quadrature>
constexpr size_t minimum_number_of_points
 Minimum number of possible collocation points for a quadrature type. More...
 
template<Basis basis>
constexpr size_t maximum_number_of_points
 Maximum number of allowed collocation points. More...
 

Detailed Description

Functionality associated with a particular choice of basis functions and quadrature for spectral operations.

Details

The functions in this namespace provide low-level access to collocation points, quadrature weights and associated matrices, such as for differentiation and interpolation. They are available in two versions: either templated directly on the enum cases of the Spectral::Basis and Spectral::Quadrature types, or taking a one-dimensional Mesh as their argument.

Note
Generally you should prefer working with a Mesh. Use its Mesh::slice_through() method to retrieve the mesh in a particular dimension:
const Mesh<2> mesh2d{
{{3, 4}},
{{Spectral::Basis::Legendre, Spectral::Basis::Legendre}},
{{Spectral::Quadrature::Gauss, Spectral::Quadrature::GaussLobatto}}};
const auto collocation_points_in_first_dim =
Spectral::collocation_points(mesh2d.slice_through(0));
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:53
const DataVector & collocation_points(size_t num_points)
Collocation points.

Most algorithms in this namespace are adapted from [109].

Enumeration Type Documentation

◆ Basis

enum class Spectral::Basis : uint8_t
strong

Either the basis functions used by a spectral or discontinuous Galerkin (DG) method, or the value FiniteDifference when a finite difference method is used.

Details

The particular choices of Basis and Quadrature determine where the collocation points of a Mesh are located in an Element. For a spectral or DG method, the Basis also represents the choice of basis functions used to represent a function on an Element, which then provides a convenient choice for the operators used for differentiation, interpolation, etc. For a finite difference method, one needs to choose the order of the scheme (and hence the weights, differentiation matrix, integration weights, and interpolant) locally in space and time to handle discontinuous solutions.

Note
Choose Legendre for a general-purpose spectral or DG mesh, unless you have a particular reason for choosing Chebyshev.
Choose two consecutive dimensions to have SphericalHarmonic to choose a spherical harmonic basis. By convention, the first dimension represents the polar/zentith angle (or colatitude), while the second dimension represents the azimuthal angle (or longitude)
We store these effectively as a 4-bit integer using the highest 4 bits of a uint8_t, which is why we do the left shift. We cannot have more than 16 bases to fit into the 4 bits, including the Uninitialized value. The number of bits to shift is encoded in the variable Spectral::detail::basis_shift.

◆ Quadrature

enum class Spectral::Quadrature : uint8_t
strong

Either the choice of quadrature method to compute integration weights for a spectral or discontinuous Galerkin (DG) method, or the locations of grid points when a finite difference method is used.

Details

The particular choices of Basis and Quadrature determine where the collocation points of a Mesh are located in an Element. For a spectral or DG method, integrals using \(N\) collocation points with Gauss quadrature are exact to polynomial order \(p=2N-1\). Gauss-Lobatto quadrature is exact only to polynomial order \(p=2N-3\), but includes collocation points at the domain boundary. For a finite difference method, one needs to choose the order of the scheme (and hence the weights, differentiation matrix, integration weights, and interpolant) locally in space and time to handle discontinuous solutions.

Note
Choose Gauss or GaussLobatto when using Basis::Legendre or Basis::Chebyshev.
Choose CellCentered or FaceCentered when using Basis::FiniteDifference.
When using Basis::SphericalHarmonic in consecutive dimensions, choose Gauss for the first dimension and Equiangular in the second dimension.
We store these effectively as a 4-bit integer using the lowest 4 bits of a uint8_t. Unlike Basis, this does not need a bitshift. We cannot have more than 16 quadratures to fit into the 4 bits, including the Uninitialized value.

Function Documentation

◆ boundary_interpolation_matrices() [1/2]

const std::pair< Matrix, Matrix > & Spectral::boundary_interpolation_matrices ( const Mesh< 1 > &  mesh)

Matrices that interpolate to the lower and upper boundaries of the element.

Assumes that the logical coordinates are \([-1, 1]\). The first element of the pair interpolates to \(\xi=-1\) and the second to \(\xi=1\). These are just the Lagrange interpolating polynomials evaluated at \(\xi=\pm1\). For Gauss-Lobatto points the only non-zero element is at the boundaries and is one and so is not implemented.

Warning
This can only be called with Gauss points.

◆ boundary_interpolation_matrices() [2/2]

template<Basis BasisType, Quadrature QuadratureType>
const std::pair< Matrix, Matrix > & Spectral::boundary_interpolation_matrices ( size_t  num_points)

Matrices that interpolate to the lower and upper boundaries of the element.

Assumes that the logical coordinates are \([-1, 1]\). The first element of the pair interpolates to \(\xi=-1\) and the second to \(\xi=1\). These are just the Lagrange interpolating polynomials evaluated at \(\xi=\pm1\). For Gauss-Lobatto points the only non-zero element is at the boundaries and is one and so is not implemented.

Warning
This can only be called with Gauss points.

◆ boundary_interpolation_term() [1/2]

const std::pair< DataVector, DataVector > & Spectral::boundary_interpolation_term ( const Mesh< 1 > &  mesh)

Interpolates values from the boundary into the volume, which is needed when applying time derivative or Bjorhus-type boundary conditions in a discontinuous Galerkin scheme using Gauss points.

Assumes that the logical coordinates are \([-1, 1]\). The interpolation is done by assuming the time derivative correction is zero on interior nodes. With a nodal Lagrange polynomial basis this means that only the \(\ell^{\mathrm{Gauss-Lobatto}}_{0}\) and \(\ell^{\mathrm{Gauss-Lobatto}}_{N}\) polynomials/basis functions contribute to the correction. In order to interpolate the correction from the nodes at the boundary, the Gauss-Lobatto Lagrange polynomials must be evaluated at the Gauss grid points. The returned pair of DataVectors stores

\begin{align*} &\ell^{\mathrm{Gauss-Lobatto}}_{0}(\xi_j^{\mathrm{Gauss}}), \\ &\ell^{\mathrm{Gauss-Lobatto}}_{N}(\xi_j^{\mathrm{Gauss}}). \end{align*}

This is a different correction from lifting. Lifting is done using the mass matrix, which is an integral over the basis functions, while here we use interpolation.

Warning
This can only be called with Gauss points.

◆ boundary_interpolation_term() [2/2]

template<Basis BasisType, Quadrature QuadratureType>
const std::pair< DataVector, DataVector > & Spectral::boundary_interpolation_term ( size_t  num_points)

Interpolates values from the boundary into the volume, which is needed when applying time derivative or Bjorhus-type boundary conditions in a discontinuous Galerkin scheme using Gauss points.

Assumes that the logical coordinates are \([-1, 1]\). The interpolation is done by assuming the time derivative correction is zero on interior nodes. With a nodal Lagrange polynomial basis this means that only the \(\ell^{\mathrm{Gauss-Lobatto}}_{0}\) and \(\ell^{\mathrm{Gauss-Lobatto}}_{N}\) polynomials/basis functions contribute to the correction. In order to interpolate the correction from the nodes at the boundary, the Gauss-Lobatto Lagrange polynomials must be evaluated at the Gauss grid points. The returned pair of DataVectors stores

\begin{align*} &\ell^{\mathrm{Gauss-Lobatto}}_{0}(\xi_j^{\mathrm{Gauss}}), \\ &\ell^{\mathrm{Gauss-Lobatto}}_{N}(\xi_j^{\mathrm{Gauss}}). \end{align*}

This is a different correction from lifting. Lifting is done using the mass matrix, which is an integral over the basis functions, while here we use interpolation.

Warning
This can only be called with Gauss points.

◆ boundary_lifting_term() [1/2]

const std::pair< DataVector, DataVector > & Spectral::boundary_lifting_term ( const Mesh< 1 > &  mesh)

Terms used during the lifting portion of a discontinuous Galerkin scheme when using Gauss points.

Assumes that the logical coordinates are \([-1, 1]\). The first element of the pair is the Lagrange polyonmials evaluated at \(\xi=-1\) divided by the weights and the second at \(\xi=1\). Specifically,

\begin{align*} \frac{\ell_j(\xi=\pm1)}{w_j} \end{align*}

Warning
This can only be called with Gauss points.

◆ boundary_lifting_term() [2/2]

template<Basis BasisType, Quadrature QuadratureType>
const std::pair< DataVector, DataVector > & Spectral::boundary_lifting_term ( size_t  num_points)

Terms used during the lifting portion of a discontinuous Galerkin scheme when using Gauss points.

Assumes that the logical coordinates are \([-1, 1]\). The first element of the pair is the Lagrange polyonmials evaluated at \(\xi=-1\) divided by the weights and the second at \(\xi=1\). Specifically,

\begin{align*} \frac{\ell_j(\xi=\pm1)}{w_j} \end{align*}

Warning
This can only be called with Gauss points.

◆ collocation_points() [1/2]

const DataVector & Spectral::collocation_points ( const Mesh< 1 > &  mesh)

Collocation points for a one-dimensional mesh.

See also
collocation_points(size_t)

◆ collocation_points() [2/2]

template<Basis BasisType, Quadrature QuadratureType>
const DataVector & Spectral::collocation_points ( size_t  num_points)

Collocation points.

Parameters
num_pointsThe number of collocation points

◆ compute_collocation_points_and_weights()

template<Basis BasisType, Quadrature QuadratureType>
std::pair< DataVector, DataVector > Spectral::compute_collocation_points_and_weights ( size_t  num_points)

Compute the collocation points and weights associated to the basis and quadrature.

Details

This function is expected to return the tuple \((\xi_k,w_k)\) where the \(\xi_k\) are the collocation points and the \(w_k\) are defined in the description of quadrature_weights(size_t).

Warning
for a FiniteDifference basis or CellCentered and FaceCentered quadratures, the weights are defined to integrate with the midpoint method

◆ differentiation_matrix() [1/2]

const Matrix & Spectral::differentiation_matrix ( const Mesh< 1 > &  mesh)

Differentiation matrix for a one-dimensional mesh.

See also
differentiation_matrix(size_t)

◆ differentiation_matrix() [2/2]

template<Basis BasisType, Quadrature QuadratureType>
const Matrix & Spectral::differentiation_matrix ( size_t  num_points)

Matrix used to compute the derivative of a function.

Details

For a function represented by the nodal coefficients \(u_j\) a matrix multiplication with the differentiation matrix \(D_{ij}\) gives the coefficients of the function's derivative. Since \(u(x)\) is expanded in Lagrange polynomials \(u(x)=\sum_j u_j l_j(x)\) the differentiation matrix is computed as \(D_{ij}=l_j^\prime(\xi_i)\) where the \(\xi_i\) are the collocation points.

The finite difference matrix uses summation by parts operators, \(D_{2-1}, D_{4-2}, D_{4-3}\), and \(D_{6-5}\) from [53].

Parameters
num_pointsThe number of collocation points

◆ differentiation_matrix_transpose() [1/2]

const Matrix & Spectral::differentiation_matrix_transpose ( const Mesh< 1 > &  mesh)

Differentiation matrix for a one-dimensional mesh.

See also
differentiation_matrix(size_t)

◆ differentiation_matrix_transpose() [2/2]

template<Basis BasisType, Quadrature QuadratureType>
const Matrix & Spectral::differentiation_matrix_transpose ( size_t  num_points)

Matrix used to compute the derivative of a function.

Details

For a function represented by the nodal coefficients \(u_j\) a matrix multiplication with the differentiation matrix \(D_{ij}\) gives the coefficients of the function's derivative. Since \(u(x)\) is expanded in Lagrange polynomials \(u(x)=\sum_j u_j l_j(x)\) the differentiation matrix is computed as \(D_{ij}=l_j^\prime(\xi_i)\) where the \(\xi_i\) are the collocation points.

The finite difference matrix uses summation by parts operators, \(D_{2-1}, D_{4-2}, D_{4-3}\), and \(D_{6-5}\) from [53].

Parameters
num_pointsThe number of collocation points

◆ integration_matrix() [1/2]

const Matrix & Spectral::integration_matrix ( const Mesh< 1 > &  mesh)

Indefinite integration matrix for a one-dimensional mesh.

See also
integration_matrix(size_t)

◆ integration_matrix() [2/2]

template<Basis BasisType, Quadrature QuadratureType>
const Matrix & Spectral::integration_matrix ( size_t  num_points)

Matrix used to perform an indefinite integral of a function over the logical grid. The left boundary condition is such that the integral is 0 at \(\xi=-1\)

Currently only Legendre and Chebyshev polynomials are implemented, but we provide a derivation for how to compute the indefinite integration matrix for general Jacobi polynomials.

Legendre Polynomials

The Legendre polynomials have the identity:

\begin{align*} P_n(x) = \frac{1}{2n+1}\frac{d}{dx}\left(P_{n+1}(x) - P_{n-1}(x)\right) \end{align*}

The goal is to evaluate the integral of a function \(u\) expanded in terms of Legendre polynomials as

\begin{align*} u(x) = \sum_{i=0}^{N} c_i P_i(x) \end{align*}

We similarly expand the indefinite integral of \(u\) as

\begin{align*} \left.\int u(y) dy\right\rvert_{x}=&\sum_{i=0}^N \tilde{c}_i P_i(x) \\ =&\left.\int\sum_{i=1}^{N}\frac{c_i}{2i+1} \left(P_{i+1}(y)-P_{i-1}(y)\right)dy\right\rvert_{x} + \tilde{c}_0 P_0(x) \\ =&\sum_{i=1}^{N}\left(\frac{c_{i-1}}{2i-1} - \frac{c_{i+1}}{2i+3}\right) P_i(x) + \tilde{c}_0 P_0(x) \end{align*}

Thus we get that for \(i>0\)

\begin{align*} \tilde{c}_i=\frac{c_{i-1}}{2i-1}-\frac{c_{i+1}}{2i+3} \end{align*}

and \(\tilde{c}_0\) is a constant of integration, which we choose such that the integral is 0 at the left boundary of the domain ( \(x=-1\)). The condition for this is:

\begin{align*} \tilde{c}_0=\sum_{i=1}^{N}(-1)^{i+1}\tilde{c}_i \end{align*}

The matrix returned by this function is the product of the tridiagonal matrix for the \(\tilde{c}_i\) and the matrix for the boundary condition.

Chebyshev Polynomials

A similar derivation leads to the relations:

\begin{align*} \tilde{c}_i=&\frac{c_{i-1}-c_{i+1}}{2i},&\mathrm{if}\;i>1 \\ \tilde{c}_1=&c_0 - \frac{c_2}{2},&\mathrm{if}\;i=1 \\ \end{align*}

We again have:

\begin{align*} \tilde{c}_0=\sum_{i=1}^N(-1)^{i+1}\tilde{c}_i \end{align*}

These are then used to define the indefinite integration matrix.

Jacobi Polynomials

For general Jacobi polynomials \(P^{(\alpha,\beta)}_n(x)\) given by

\begin{align*} (1-x)^\alpha(1+x)^\beta P^{(\alpha,\beta)}_n(x)=\frac{(-1)^n}{2^n n!} \frac{d^n}{dx^n}\left[(1-x)^{\alpha+n}(1+x)^{\beta+n}\right] \end{align*}

we have that

\begin{align*} P^{(\alpha,\beta)}_n(x)=\frac{d}{dx}\left( b^{(\alpha,\beta)}_{n-1,n}P^{(\alpha,\beta)}_{n-1}(x) + b^{(\alpha,\beta)}_{n,n}P^{(\alpha,\beta)}_n(x) + b^{(\alpha,\beta)}_{n+1,n}P^{(\alpha,\beta)}_{n+1}(x) \right) \end{align*}

where

\begin{align*} b^{(\alpha,\beta)}_{n-1,n}=&-\frac{1}{n+\alpha+\beta} a^{(\alpha,\beta)}_{n-1,n} \\ b^{(\alpha,\beta)}_{n,n}=&-\frac{2}{\alpha+\beta} a^{(\alpha,\beta)}_{n,n} \\ b^{(\alpha,\beta)}_{n+1,n}=&\frac{1}{n+1} a^{(\alpha,\beta)}_{n+1,n} \\ a^{(\alpha,\beta)}_{n-1,n}=&\frac{2(n+\alpha)(n+\beta)} {(2n+\alpha+\beta+1)(2n+\alpha+\beta)} \\ a^{(\alpha,\beta)}_{n,n}=&-\frac{\alpha^2-\beta^2} {(2n+\alpha+\beta+2)(2n+\alpha+\beta)} \\ a^{(\alpha,\beta)}_{n-1,n}=&\frac{2(n+1)(n+\alpha+\beta+1)} {(2n+\alpha+\beta+2)(2n+\alpha+\beta+1)} \end{align*}

Following the same derivation we get that

\begin{align*} \tilde{c}_i=c_{i+1}b^{(\alpha,\beta)}_{i,i+1} +c_i b^{(\alpha,\beta)}_{i,i} +c_{i-1}b^{(\alpha,\beta)}_{i,i-1} \end{align*}

and the boundary condition is

\begin{align*} \tilde{c}_0=\sum_{i=1}^N(-1)^{i+1} \frac{\Gamma(i+\alpha+1)}{i!\Gamma(\alpha+1)} \tilde{c}_i \end{align*}

where \(\Gamma(x)\) is the Gamma function.

◆ interpolation_matrix() [1/2]

template<typename T >
Matrix Spectral::interpolation_matrix ( const Mesh< 1 > &  mesh,
const T &  target_points 
)

Interpolation matrix to the target_points for a one-dimensional mesh.

See also
interpolation_matrix(size_t, const T&)

◆ interpolation_matrix() [2/2]

template<Basis BasisType, Quadrature QuadratureType, typename T >
Matrix Spectral::interpolation_matrix ( size_t  num_points,
const T &  target_points 
)

Matrix used to interpolate to the target_points.

Warning
For each target point located outside of the logical coordinate interval covered by BasisType (often \([-1,1]\)), the resulting matrix performs polynomial extrapolation instead of interpolation. The extrapolation will be correct but may suffer from reduced accuracy, especially for higher-order polynomials (i.e., larger values of num_points).
Parameters
num_pointsThe number of collocation points
target_pointsThe points to interpolate to

◆ linear_filter_matrix() [1/2]

const Matrix & Spectral::linear_filter_matrix ( const Mesh< 1 > &  mesh)

Linear filter matrix for a one-dimensional mesh.

See also
linear_filter_matrix(size_t)

◆ linear_filter_matrix() [2/2]

template<Basis BasisType, Quadrature QuadratureType>
const Matrix & Spectral::linear_filter_matrix ( size_t  num_points)

Matrix used to linearize a function.

Details

Filters out all except the lowest two modes by applying \(\mathcal{V}^{-1}\cdot\mathrm{diag}(1,1,0,0,...)\cdot\mathcal{V}\) to the nodal coefficients, where \(\mathcal{V}\) is the Vandermonde matrix computed in modal_to_nodal_matrix(size_t).

Parameters
num_pointsThe number of collocation points
See also
modal_to_nodal_matrix(size_t)
nodal_to_modal_matrix(size_t)

◆ modal_to_nodal_matrix() [1/2]

const Matrix & Spectral::modal_to_nodal_matrix ( const Mesh< 1 > &  mesh)

Transformation matrix from modal to nodal coefficients for a one-dimensional mesh.

See also
modal_to_nodal_matrix(size_t)

◆ modal_to_nodal_matrix() [2/2]

template<Basis BasisType, Quadrature QuadratureType>
const Matrix & Spectral::modal_to_nodal_matrix ( size_t  num_points)

Matrix used to transform from the spectral coefficients (modes) of a function to its nodal coefficients. Also referred to as the Vandermonde matrix.

Details

The Vandermonde matrix is computed as \(\mathcal{V}_{ij}=\Phi_j(\xi_i)\) where the \(\Phi_j(x)\) are the spectral basis functions used in the modal expansion \(u(x)=\sum_j \widetilde{u}_j\Phi_j(x)\), e.g. normalized Legendre polynomials, and the \(\xi_i\) are the collocation points used to construct the interpolating Lagrange polynomials in the nodal expansion \(u(x)=\sum_j u_j l_j(x)\). Then the Vandermonde matrix arises as \(u(\xi_i)=u_i=\sum_j \widetilde{u}_j\Phi_j(\xi_i)=\sum_j \mathcal{V}_{ij}\widetilde{u}_j\).

Parameters
num_pointsThe number of collocation points
See also
nodal_to_modal_matrix(size_t)

◆ nodal_to_modal_matrix() [1/2]

const Matrix & Spectral::nodal_to_modal_matrix ( const Mesh< 1 > &  mesh)

Transformation matrix from nodal to modal coefficients for a one-dimensional mesh.

See also
nodal_to_modal_matrix(size_t)

◆ nodal_to_modal_matrix() [2/2]

template<Basis BasisType, Quadrature QuadratureType>
const Matrix & Spectral::nodal_to_modal_matrix ( size_t  num_points)

Matrix used to transform from the nodal coefficients of a function to its spectral coefficients (modes). Also referred to as the inverse Vandermonde matrix.

Details

This is the inverse to the Vandermonde matrix \(\mathcal{V}\) computed in modal_to_nodal_matrix(size_t). It can be computed analytically for Gauss quadrature by evaluating \(\sum_j\mathcal{V}^{-1}_{ij}u_j=\widetilde{u}_i= \frac{(u,\Phi_i)}{\gamma_i}\) for a Lagrange basis function \(u(x)=l_k(x)\) to find \(\mathcal{V}^{-1}_{ij}=\mathcal{V}_{ji}\frac{w_j}{\gamma_i}\) where the \(w_j\) are the Gauss quadrature weights and \(\gamma_i\) is the norm square of the spectral basis function \(\Phi_i\).

Parameters
num_pointsThe number of collocation points
See also
modal_to_nodal_matrix(size_t)

◆ projection_matrix_child_to_parent()

const Matrix & Spectral::projection_matrix_child_to_parent ( const Mesh< 1 > &  child_mesh,
const Mesh< 1 > &  parent_mesh,
ChildSize  size,
bool  operand_is_massive = false 
)

The projection matrix from a child mesh to its parent.

The projection matrices returned by this function (and by projection_matrix_parent_to_child()) define orthogonal projection operators between the spaces of functions on a parent mesh and its children. These projections are usually the correct way to transfer data between meshes in a mesh-refinement hierarchy, as well as between an element face and its adjacent mortars.

These functions assume that the child_mesh is at least as fine as the parent_mesh, i.e. functions on the parent_mesh can be represented exactly on the child_mesh. In practice this means that functions can be projected to a mortar (the child_mesh) from both adjacent element faces (the parent_mesh) without losing accuracy. Similarly, functions in a mesh-refinement hierarchy don't lose accuracy when an element is split (h-refined). For this reason, the projection_matrix_child_to_parent is sometimes referred to as a "restriction operator" and the projection_matrix_parent_to_child as a "prolongation operator".

Massive quantities
If the quantity that should be projected is not a function over the computational grid but a "massive" residual, i.e. a quantity \(\int_{\Omega_k} f(x) \psi_p(x) \mathrm{d}V\) where \(\psi_p\) are the basis functions on the mesh, then pass true for the parameter operand_is_massive (default is false). The restriction operator for this case is just the transpose of the prolongation operator, i.e. just an interpolation matrix transpose. Note that the "massive" residual already takes the difference in element size between parent and children into account by including a Jacobian in the volume element of the integral.
Implementation details
The half-interval projections are based on an equation derived by Saul. This shows that the projection from the spectral basis for the entire interval to the spectral basis for the upper half interval is

\begin{equation*} T_{jk} = \frac{2 j + 1}{2} 2^j \sum_{n=0}^{j-k} \binom{j}{k+n} \binom{(j + k + n - 1)/2}{j} \frac{(k + n)!^2}{(2 k + n + 1)! n!} \end{equation*}

◆ projection_matrix_parent_to_child()

const Matrix & Spectral::projection_matrix_parent_to_child ( const Mesh< 1 > &  parent_mesh,
const Mesh< 1 > &  child_mesh,
ChildSize  size 
)

The projection matrix from a parent mesh to one of its children.

See also
projection_matrix_child_to_parent()

◆ quadrature_weights() [1/2]

const DataVector & Spectral::quadrature_weights ( const Mesh< 1 > &  mesh)

Quadrature weights for a one-dimensional mesh.

See also
quadrature_weights(size_t)

◆ quadrature_weights() [2/2]

template<Basis BasisType, Quadrature QuadratureType>
const DataVector & Spectral::quadrature_weights ( size_t  num_points)

Weights to compute definite integrals.

Details

These are the coefficients to contract with the nodal function values \(f_k\) to approximate the definite integral \(I[f]=\int f(x)\mathrm{d}x\).

Note that the term quadrature also often refers to the quantity \(Q[f]=\int f(x)w(x)\mathrm{d}x\approx \sum_k f_k w_k\). Here, \(w(x)\) denotes the basis-specific weight function w.r.t. to which the basis functions \(\Phi_k\) are orthogonal, i.e \(\int\Phi_i(x)\Phi_j(x)w(x)=0\) for \(i\neq j\). The weights \(w_k\) approximate this inner product. To approximate the definite integral \(I[f]\) we must employ the coefficients \(\frac{w_k}{w(\xi_k)}\) instead, where the \(\xi_k\) are the collocation points. These are the coefficients this function returns. Only for a unit weight function \(w(x)=1\), i.e. a Legendre basis, is \(I[f]=Q[f]\) so this function returns the \(w_k\) identically.

For a FiniteDifference basis or CellCentered and FaceCentered quadratures, the interpretation of the quadrature weights in term of an approximation to \(I(q)\) remains correct, but its explanation in terms of orthonormal basis is not, i.e. we set \(w_k\) to the grid spacing at each point, and the inverse weight \(\frac{1}{w(\xi_k)}=1\) to recover the midpoint method for definite integrals.

Parameters
num_pointsThe number of collocation points

◆ weak_flux_differentiation_matrix() [1/2]

const Matrix & Spectral::weak_flux_differentiation_matrix ( const Mesh< 1 > &  mesh)

Matrix used to compute the divergence of the flux in weak form.

See also
weak_flux_differentiation_matrix(size_t)

◆ weak_flux_differentiation_matrix() [2/2]

template<Basis BasisType, Quadrature QuadratureType>
const Matrix & Spectral::weak_flux_differentiation_matrix ( size_t  num_points)

Matrix used to compute the divergence of the flux in weak form.

This is the transpose of the differentiation matrix multiplied by quadrature weights that appear in DG integrals:

\begin{equation} \frac{D^T_{ij}} \frac{w_j}{w_i} \end{equation}

Parameters
num_pointsThe number of collocation points

Variable Documentation

◆ basis_shift

constexpr uint8_t Spectral::basis_shift = 4
constexpr

The amount we left-shift the basis integers by.

We do this to be able to represent the combined Basis and Quadrature as one 8-bit integer.

◆ maximum_number_of_points

template<Basis basis>
constexpr size_t Spectral::maximum_number_of_points
constexpr
Initial value:
=
basis == Basis::FiniteDifference ? 40 : 20

Maximum number of allowed collocation points.

Details

We choose a limit of 24 FD grid points because for DG-subcell the number of points in an element is 2 * number_dg_points - 1 for cell centered, and 2 * number_dg_points for face-centered. Because there is no way of generically retrieving the maximum number of grid points for a non-FD basis, we need to hard-code both values here. If the number of grid points is increased for the non-FD bases, it should also be increased for the FD basis. Note that for good task-based parallelization 24 grid points is already a fairly large number.

◆ minimum_number_of_points

template<Basis basis, Quadrature quadrature>
constexpr size_t Spectral::minimum_number_of_points
constexpr
Initial value:
=
detail::minimum_number_of_points(basis, quadrature)

Minimum number of possible collocation points for a quadrature type.

Details

Since Gauss-Lobatto quadrature has points on the domain boundaries it must have at least two collocation points. Gauss quadrature can have only one collocation point.

Details

For CellCentered the minimum number of points is 1, while for FaceCentered it is 2.