Spectral.hpp
Go to the documentation of this file.
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 /*!
5  * \file
6  * Declares functionality to retrieve spectral quantities associated with
7  * a particular choice of basis functions and quadrature.
8  */
9 
10 #pragma once
11 
12 #include <cstddef>
13 #include <iosfwd>
14 #include <limits>
15 #include <utility>
16 
17 /// \cond
18 class Matrix;
19 class DataVector;
20 template <size_t>
21 class Mesh;
22 /// \endcond
23 
24 /*!
25  * \ingroup SpectralGroup
26  * \brief Functionality associated with a particular choice of basis functions
27  * and quadrature for spectral operations.
28  *
29  * \details The functions in this namespace provide low-level access to
30  * collocation points, quadrature weights and associated matrices, such as for
31  * differentiation and interpolation. They are available in two versions: either
32  * templated directly on the enum cases of the Spectral::Basis and
33  * Spectral::Quadrature types, or taking a one-dimensional Mesh as their
34  * argument.
35  *
36  * \note Generally you should prefer working with a Mesh. Use its
37  * Mesh::slice_through() method to retrieve the mesh in a particular dimension:
38  * \snippet Test_Spectral.cpp get_points_for_mesh
39  *
40  *
41  * Most algorithms in this namespace are adapted from \cite Kopriva.
42  */
43 namespace Spectral {
44 
45 /*!
46  * \ingroup SpectralGroup
47  * \brief The choice of basis functions for computing collocation points and
48  * weights.
49  *
50  * \details Choose `Legendre` for a general-purpose DG mesh, unless you have a
51  * particular reason for choosing another basis.
52  */
53 enum class Basis { Chebyshev, Legendre };
54 
55 /// \cond HIDDEN_SYMBOLS
56 std::ostream& operator<<(std::ostream& os, const Basis& basis) noexcept;
57 /// \endcond
58 
59 /*!
60  * \ingroup SpectralGroup
61  * \brief The choice of quadrature method to compute integration weights.
62  *
63  * \details Integrals using \f$N\f$ collocation points with Gauss quadrature are
64  * exact to polynomial order \f$p=2N-1\f$. Gauss-Lobatto quadrature is exact
65  * only to polynomial order \f$p=2N-3\f$, but includes collocation points at the
66  * domain boundary.
67  */
68 enum class Quadrature { Gauss, GaussLobatto };
69 
70 /// \cond HIDDEN_SYMBOLS
71 std::ostream& operator<<(std::ostream& os,
72  const Quadrature& quadrature) noexcept;
73 /// \endcond
74 
75 /*!
76  * \brief Minimum number of possible collocation points for a quadrature type.
77  *
78  * \details Since Gauss-Lobatto quadrature has points on the domain boundaries
79  * it must have at least two collocation points. Gauss quadrature can have only
80  * one collocation point.
81  */
82 template <Basis, Quadrature>
84 /// \cond
85 template <Basis BasisType>
86 constexpr size_t minimum_number_of_points<BasisType, Quadrature::Gauss> = 1;
87 template <Basis BasisType>
88 constexpr size_t minimum_number_of_points<BasisType, Quadrature::GaussLobatto> =
89  2;
90 /// \endcond
91 
92 /*!
93  * \brief Maximum number of allowed collocation points.
94  */
95 template <Basis>
96 constexpr size_t maximum_number_of_points = 12;
97 
98 /*!
99  * \brief Compute the function values of the basis function \f$\Phi_k(x)\f$
100  * (zero-indexed).
101  */
102 template <Basis BasisType, typename T>
103 T compute_basis_function_value(size_t k, const T& x) noexcept;
104 
105 /*!
106  * \brief Compute the inverse of the weight function \f$w(x)\f$ w.r.t. which
107  * the basis functions are orthogonal. See the description of
108  * `quadrature_weights(size_t)` for details.
109  */
110 template <Basis>
112 
113 /*!
114  * \brief Compute the normalization square of the basis function \f$\Phi_k\f$
115  * (zero-indexed), i.e. the weighted definite integral over its square.
116  */
117 template <Basis BasisType>
118 double compute_basis_function_normalization_square(size_t k) noexcept;
119 
120 /*!
121  * \brief Compute the collocation points and weights associated to the
122  * basis and quadrature.
123  *
124  * \details This function is expected to return the tuple
125  * \f$(\xi_k,w_k)\f$ where the \f$\xi_k\f$ are the collocation
126  * points and the \f$w_k\f$ are defined in the description of
127  * `quadrature_weights(size_t)`.
128  */
129 template <Basis BasisType, Quadrature QuadratureType>
131  size_t num_points) noexcept;
132 
133 /*!
134  * \brief Collocation points
135  * \param num_points The number of collocation points
136  */
137 template <Basis BasisType, Quadrature QuadratureType>
138 const DataVector& collocation_points(size_t num_points) noexcept;
139 
140 /*!
141  * \brief Collocation points for a one-dimensional mesh.
142  *
143  * \see collocation_points(size_t)
144  */
145 const DataVector& collocation_points(const Mesh<1>& mesh) noexcept;
146 
147 /*!
148  * \brief Weights to compute definite integrals.
149  *
150  * \details These are the coefficients to contract with the nodal
151  * function values \f$f_k\f$ to approximate the definite integral \f$I[f]=\int
152  * f(x)\mathrm{d}x\f$.
153  *
154  * Note that the term _quadrature_ also often refers to the quantity
155  * \f$Q[f]=\int f(x)w(x)\mathrm{d}x\approx \sum_k f_k w_k\f$. Here, \f$w(x)\f$
156  * denotes the basis-specific weight function w.r.t. to which the basis
157  * functions \f$\Phi_k\f$ are orthogonal, i.e \f$\int\Phi_i(x)\Phi_j(x)w(x)=0\f$
158  * for \f$i\neq j\f$. The weights \f$w_k\f$ approximate this inner product. To
159  * approximate the definite integral \f$I[f]\f$ we must employ the
160  * coefficients \f$\frac{w_k}{w(\xi_k)}\f$ instead, where the \f$\xi_k\f$ are
161  * the collocation points. These are the coefficients this function returns.
162  * Only for a unit weight function \f$w(x)=1\f$, i.e. a Legendre basis, is
163  * \f$I[f]=Q[f]\f$ so this function returns the \f$w_k\f$ identically.
164  *
165  * \param num_points The number of collocation points
166  */
167 template <Basis BasisType, Quadrature QuadratureType>
168 const DataVector& quadrature_weights(size_t num_points) noexcept;
169 
170 /*!
171  * \brief Quadrature weights for a one-dimensional mesh.
172  *
173  * \see quadrature_weights(size_t)
174  */
175 const DataVector& quadrature_weights(const Mesh<1>& mesh) noexcept;
176 
177 /*!
178  * \brief %Matrix used to compute the derivative of a function.
179  *
180  * \details For a function represented by the nodal coefficients \f$u_j\f$ a
181  * matrix multiplication with the differentiation matrix \f$D_{ij}\f$ gives the
182  * coefficients of the function's derivative. Since \f$u(x)\f$ is expanded in
183  * Lagrange polynomials \f$u(x)=\sum_j u_j l_j(x)\f$ the differentiation matrix
184  * is computed as \f$D_{ij}=l_j^\prime(\xi_i)\f$ where the \f$\xi_i\f$ are the
185  * collocation points.
186  *
187  * \param num_points The number of collocation points
188  */
189 template <Basis BasisType, Quadrature QuadratureType>
190 const Matrix& differentiation_matrix(size_t num_points) noexcept;
191 
192 /*!
193  * \brief Differentiation matrix for a one-dimensional mesh.
194  *
195  * \see differentiation_matrix(size_t)
196  */
197 const Matrix& differentiation_matrix(const Mesh<1>& mesh) noexcept;
198 
199 /*!
200  * \brief %Matrix used to interpolate to the \p target_points.
201  *
202  * \warning For each target point located outside of the logical coordinate
203  * interval covered by `BasisType` (often \f$[-1,1]\f$), the resulting matrix
204  * performs polynomial extrapolation instead of interpolation. The extapolation
205  * will be correct but may suffer from reduced accuracy, especially for
206  * higher-order polynomials (i.e., larger values of `num_points`).
207  *
208  * \param num_points The number of collocation points
209  * \param target_points The points to interpolate to
210  */
211 template <Basis BasisType, Quadrature QuadratureType, typename T>
212 Matrix interpolation_matrix(size_t num_points, const T& target_points) noexcept;
213 
214 /*!
215  * \brief Interpolation matrix to the \p target_points for a one-dimensional
216  * mesh.
217  *
218  * \see interpolation_matrix(size_t, const T&)
219  */
220 template <typename T>
222  const T& target_points) noexcept;
223 
224 /*!
225  * \brief %Matrix used to transform from the spectral coefficients (modes) of a
226  * function to its nodal coefficients. Also referred to as the _Vandermonde
227  * matrix_.
228  *
229  * \details The Vandermonde matrix is computed as
230  * \f$\mathcal{V}_{ij}=\Phi_j(\xi_i)\f$ where the \f$\Phi_j(x)\f$ are the
231  * spectral basis functions used in the modal expansion
232  * \f$u(x)=\sum_j \widetilde{u}_j\Phi_j(x)\f$, e.g. normalized Legendre
233  * polynomials, and the \f$\xi_i\f$ are the collocation points used to construct
234  * the interpolating Lagrange polynomials in the nodal expansion
235  * \f$u(x)=\sum_j u_j l_j(x)\f$. Then the Vandermonde matrix arises as
236  * \f$u(\xi_i)=u_i=\sum_j \widetilde{u}_j\Phi_j(\xi_i)=\sum_j
237  * \mathcal{V}_{ij}\widetilde{u}_j\f$.
238  *
239  * \param num_points The number of collocation points
240 
241  * \see grid_points_to_spectral_matrix(size_t)
242  */
243 template <Basis BasisType, Quadrature QuadratureType>
244 const Matrix& spectral_to_grid_points_matrix(size_t num_points) noexcept;
245 
246 /*!
247  * \brief Transformation matrix from modal to nodal coefficients for a
248  * one-dimensional mesh.
249  *
250  * \see spectral_to_grid_points_matrix(size_t)
251  */
252 const Matrix& spectral_to_grid_points_matrix(const Mesh<1>& mesh) noexcept;
253 
254 /*!
255  * \brief %Matrix used to transform from the nodal coefficients of a function to
256  * its spectral coefficients (modes). Also referred to as the inverse
257  * _Vandermonde matrix_.
258  *
259  * \details This is the inverse to the Vandermonde matrix \f$\mathcal{V}\f$
260  * computed in spectral_to_grid_points_matrix(size_t). It can be computed
261  * analytically for Gauss quadrature by evaluating
262  * \f$\sum_j\mathcal{V}^{-1}_{ij}u_j=\widetilde{u}_i=
263  * \frac{(u,\Phi_i)}{\gamma_i}\f$
264  * for a Lagrange basis function \f$u(x)=l_k(x)\f$ to find
265  * \f$\mathcal{V}^{-1}_{ij}=\mathcal{V}_{ji}\frac{w_j}{\gamma_i}\f$ where the
266  * \f$w_j\f$ are the Gauss quadrature weights and \f$\gamma_i\f$ is the norm
267  * square of the spectral basis function \f$\Phi_i\f$.
268  *
269  * \param num_points The number of collocation points
270  *
271  * \see spectral_to_grid_points_matrix(size_t)
272  */
273 template <Basis BasisType, Quadrature QuadratureType>
274 const Matrix& grid_points_to_spectral_matrix(size_t num_points) noexcept;
275 
276 /*!
277  * \brief Transformation matrix from nodal to modal coefficients for a
278  * one-dimensional mesh.
279  *
280  * \see grid_points_to_spectral_matrix(size_t)
281  */
282 const Matrix& grid_points_to_spectral_matrix(const Mesh<1>& mesh) noexcept;
283 
284 /*!
285  * \brief %Matrix used to linearize a function.
286  *
287  * \details Filters out all except the lowest two modes by applying
288  * \f$\mathcal{V}^{-1}\cdot\mathrm{diag}(1,1,0,0,...)\cdot\mathcal{V}\f$ to the
289  * nodal coefficients, where \f$\mathcal{V}\f$ is the Vandermonde matrix
290  * computed in `spectral_to_grid_points_matrix(size_t)`.
291  *
292  * \param num_points The number of collocation points
293  *
294  * \see spectral_to_grid_points_matrix(size_t)
295  * \see grid_points_to_spectral_matrix(size_t)
296  */
297 template <Basis BasisType, Quadrature QuadratureType>
298 const Matrix& linear_filter_matrix(size_t num_points) noexcept;
299 
300 /*!
301  * \brief Linear filter matrix for a one-dimensional mesh.
302  *
303  * \see linear_filter_matrix(size_t)
304  */
305 const Matrix& linear_filter_matrix(const Mesh<1>& mesh) noexcept;
306 
307 } // namespace Spectral
Quadrature
The choice of quadrature method to compute integration weights.
Definition: Spectral.hpp:68
const Matrix & spectral_to_grid_points_matrix(size_t num_points) noexcept
Matrix used to transform from the spectral coefficients (modes) of a function to its nodal coefficien...
DataVector compute_inverse_weight_function_values(const DataVector &) noexcept
Compute the inverse of the weight function w.r.t. which the basis functions are orthogonal. See the description of quadrature_weights(size_t) for details.
const DataVector & collocation_points(const size_t num_points) noexcept
Collocation points.
Definition: Spectral.cpp:232
Holds the number of grid points, basis, and quadrature in each direction of the computational grid...
Definition: Mesh.hpp:49
constexpr size_t maximum_number_of_points
Maximum number of allowed collocation points.
Definition: Spectral.hpp:96
Basis
The choice of basis functions for computing collocation points and weights.
Definition: Spectral.hpp:53
double compute_basis_function_normalization_square(size_t k) noexcept
Compute the normalization square of the basis function (zero-indexed), i.e. the weighted definite in...
Matrix interpolation_matrix(const size_t num_points, const T &target_points) noexcept
Matrix used to interpolate to the target_points.
Definition: Spectral.cpp:273
T compute_basis_function_value(size_t k, const T &x) noexcept
Compute the function values of the basis function (zero-indexed).
const Matrix & linear_filter_matrix(size_t num_points) noexcept
Matrix used to linearize a function.
const Matrix & differentiation_matrix(size_t num_points) noexcept
Matrix used to compute the derivative of a function.
T max(T... args)
A dynamically sized matrix of doubles with column-major storage.
Definition: Matrix.hpp:19
std::pair< DataVector, DataVector > compute_collocation_points_and_weights(size_t num_points) noexcept
Compute the collocation points and weights associated to the basis and quadrature.
constexpr size_t minimum_number_of_points
Minimum number of possible collocation points for a quadrature type.
Definition: Spectral.hpp:83
Stores a collection of function values.
Definition: DataVector.hpp:46
Functionality associated with a particular choice of basis functions and quadrature for spectral oper...
Definition: Chebyshev.cpp:16
const Matrix & grid_points_to_spectral_matrix(size_t num_points) noexcept
Matrix used to transform from the nodal coefficients of a function to its spectral coefficients (mode...
const DataVector & quadrature_weights(const size_t num_points) noexcept
Weights to compute definite integrals.
Definition: Spectral.cpp:241