|
SpECTRE
v2024.05.11
|
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 { Chebyshev , Legendre , FiniteDifference , SphericalHarmonic } |
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 { Full , UpperHalf , LowerHalf } |
| The portion of a mesh covered by a child mesh. | |
| enum class | Quadrature { 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, 4 > | all_bases () |
| All possible values of Basis. | |
| Basis | to_basis (const std::string &basis) |
| Convert a string to a Basis enum. | |
| std::ostream & | operator<< (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::ostream & | operator<< (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 Matrix & | 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. 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 Matrix & | 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. 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, 5 > | all_quadratures () |
| All possible values of Quadrature. | |
| Quadrature | to_quadrature (const std::string &quadrature) |
| Convert a string to a Quadrature enum. | |
| std::ostream & | operator<< (std::ostream &os, const Quadrature &quadrature) |
| Output operator for a Quadrature. | |
| template<Basis BasisType, typename T > | |
| 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, DataVector > | compute_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 DataVector & | collocation_points (size_t num_points) |
| Collocation points. More... | |
| const DataVector & | collocation_points (const Mesh< 1 > &mesh) |
| Collocation points for a one-dimensional mesh. More... | |
| template<Basis BasisType, Quadrature QuadratureType> | |
| const DataVector & | quadrature_weights (size_t num_points) |
| Weights to compute definite integrals. More... | |
| const DataVector & | quadrature_weights (const Mesh< 1 > &mesh) |
| Quadrature weights for a one-dimensional mesh. More... | |
| template<Basis BasisType, Quadrature QuadratureType> | |
| const Matrix & | weak_flux_differentiation_matrix (size_t num_points) |
| Matrix used to compute the divergence of the flux in weak form. More... | |
| const Matrix & | weak_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 Matrix & | 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\) More... | |
| const Matrix & | integration_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 Matrix & | 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. More... | |
| const Matrix & | modal_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 Matrix & | 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. More... | |
| const Matrix & | nodal_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 Matrix & | linear_filter_matrix (size_t num_points) |
| Matrix used to linearize a function. More... | |
| const Matrix & | linear_filter_matrix (const Mesh< 1 > &mesh) |
| Linear filter matrix for a one-dimensional mesh. More... | |
| template<Basis BasisType, Quadrature QuadratureType> | |
| const Matrix & | differentiation_matrix (size_t num_points) |
| Matrix used to compute the derivative of a function. More... | |
| template<Basis BasisType, Quadrature QuadratureType> | |
| const Matrix & | differentiation_matrix_transpose (size_t num_points) |
| Matrix used to compute the derivative of a function. More... | |
| const Matrix & | differentiation_matrix (const Mesh< 1 > &mesh) |
| Differentiation matrix for a one-dimensional mesh. More... | |
| const Matrix & | differentiation_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 | |
| 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... | |
Functionality associated with a particular choice of basis functions and quadrature for spectral operations.
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.
Most algorithms in this namespace are adapted from [109].
|
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.
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.
Legendre for a general-purpose spectral or DG mesh, unless you have a particular reason for choosing Chebyshev.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)
|
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.
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.
Gauss or GaussLobatto when using Basis::Legendre or Basis::Chebyshev.CellCentered or FaceCentered when using Basis::FiniteDifference.Gauss for the first dimension and Equiangular in the second dimension. | 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.
| 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.
| 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.
| 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.
| 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*}
| 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*}
| const DataVector & Spectral::collocation_points | ( | const Mesh< 1 > & | mesh | ) |
Collocation points for a one-dimensional mesh.
| const DataVector & Spectral::collocation_points | ( | size_t | num_points | ) |
Collocation points.
| num_points | The number of collocation points |
| 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.
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).
FiniteDifference basis or CellCentered and FaceCentered quadratures, the weights are defined to integrate with the midpoint method Differentiation matrix for a one-dimensional mesh.
| const Matrix & Spectral::differentiation_matrix | ( | size_t | num_points | ) |
Matrix used to compute the derivative of a function.
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].
| num_points | The number of collocation points |
Differentiation matrix for a one-dimensional mesh.
| const Matrix & Spectral::differentiation_matrix_transpose | ( | size_t | num_points | ) |
Matrix used to compute the derivative of a function.
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].
| num_points | The number of collocation points |
Indefinite integration matrix for a one-dimensional mesh.
| 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.
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.
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.
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.
| Matrix Spectral::interpolation_matrix | ( | const Mesh< 1 > & | mesh, |
| const T & | target_points | ||
| ) |
Interpolation matrix to the target_points for a one-dimensional mesh.
| Matrix Spectral::interpolation_matrix | ( | size_t | num_points, |
| const T & | target_points | ||
| ) |
Matrix used to interpolate to the target_points.
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).| num_points | The number of collocation points |
| target_points | The points to interpolate to |
Linear filter matrix for a one-dimensional mesh.
| const Matrix & Spectral::linear_filter_matrix | ( | size_t | num_points | ) |
Matrix used to linearize a function.
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).
| num_points | The number of collocation points |
Transformation matrix from modal to nodal coefficients for a one-dimensional mesh.
| 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.
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\).
| num_points | The number of collocation points |
Transformation matrix from nodal to modal coefficients for a one-dimensional mesh.
| 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.
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\).
| num_points | The number of collocation points |
| 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".
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.\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*}
| 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.
| const DataVector & Spectral::quadrature_weights | ( | const Mesh< 1 > & | mesh | ) |
Quadrature weights for a one-dimensional mesh.
| const DataVector & Spectral::quadrature_weights | ( | size_t | num_points | ) |
Weights to compute definite integrals.
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.
| num_points | The number of collocation points |
Matrix used to compute the divergence of the flux in weak form.
| 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}
| num_points | The number of collocation points |
|
constexpr |
Maximum number of allowed collocation points.
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.
|
constexpr |
Minimum number of possible collocation points for a quadrature type.
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.
For CellCentered the minimum number of points is 1, while for FaceCentered it is 2.