Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <limits>
9 :
10 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
11 : #include "Utilities/Gsl.hpp"
12 : #include "Utilities/TMPL.hpp"
13 :
14 : /// \cond
15 : class DataVector;
16 : /// \endcond
17 :
18 : /*!
19 : * \brief Items for assessing truncation error in spectral methods.
20 : */
21 1 : namespace PowerMonitors {
22 :
23 : /// @{
24 : /*!
25 : * \ingroup SpectralGroup
26 : * \brief Returns array of power monitors in each spatial dimension.
27 : *
28 : * Computed following Sec. 5.1 of Ref. \cite Szilagyi2014fna.
29 : * For example, in the x dimension (indexed by \f$ k_0 \f$), we compute
30 : *
31 : * \f{align*}{
32 : * P_{k_0}[\psi] = \sqrt{ \frac{1}{N_1 N_2}
33 : * \sum_{k_1,k_2} \left| C_{k_0,k_1,k_2} \right|^2} ,
34 : * \f}
35 : *
36 : * where \f$ C_{k_0,k_1,k_2}\f$ are the modal coefficients
37 : * of variable \f$ \psi \f$.
38 : *
39 : * For dimensions using a Fourier basis, the cosine and sine power for each
40 : * wavenumber \f$ k \f$ are combined via \f$ P_k = \sqrt{P_{\cos,k}^2 +
41 : * P_{\sin,k}^2} \f$, so the output size for a Fourier dimension of \f$ N \f$
42 : * points is \f$ N/2 + 1 \f$ (integer division) rather than \f$ N \f$.
43 : *
44 : */
45 : template <typename VectorType, size_t Dim>
46 1 : void power_monitors(gsl::not_null<std::array<DataVector, Dim>*> result,
47 : const VectorType& u, const Mesh<Dim>& mesh);
48 :
49 : template <typename VectorType, size_t Dim>
50 1 : std::array<DataVector, Dim> power_monitors(const VectorType& u,
51 : const Mesh<Dim>& mesh);
52 : /// @}
53 :
54 : /// @{
55 : /*!
56 : * \ingroup SpectralGroup
57 : * \brief Compute the relative truncation error.
58 : *
59 : * The negative logarithm of this quantity is defined by Eqs. (57) and
60 : * (58) of Ref. \cite Szilagyi2014fna, i.e.,
61 : *
62 : * \f{align*}{
63 : * \mathcal{T}\left[P_k\right] = \log_{10} \max \left(P_0, P_1\right)
64 : * - \dfrac{\sum_{j=0}^{j_{\text{max}, k}} \log_{10} \left(P_j\right) w_j}
65 : * {\sum_{j=0}^{j_{\text{max}, k}} w_j} , \f}
66 : *
67 : * with weights
68 : *
69 : * \f{align*}{
70 : * w_j = \exp\left[ - \left(j - j_{\text{max}, k}
71 : * + \dfrac{1}{2}\right)^2 \right] .
72 : * \f}
73 : *
74 : * where \f$ j_{\text{max}, k} = N_k - 1 \f$ and \f$ N_k \f$ is the number of
75 : * modes or gridpoints in dimension k. Here the second term is a weighted
76 : * average with larger weights toward the highest modes.
77 : *
78 : * \note Modes below a cutoff of $100 \epsilon \mathrm{max}_k(P_k)$ are ignored
79 : * in the weighted average, where $\epsilon$ is the machine epsilon. This
80 : * ensures that we don't underestimate the truncation error if some modes are
81 : * zero (e.g. by symmetry). Furthermore, if the last two or more modes are zero,
82 : * we assume that the function is represented exactly and return a relative
83 : * truncation error of zero.
84 : *
85 : * \details The number of modes (`num_modes_to_use`) argument needs to be less
86 : * or equal than the total number of power monitors (`power_monitor.size()`).
87 : * In contrast with Ref. \cite Szilagyi2014fna, here we index the modes starting
88 : * from zero.
89 : *
90 : */
91 1 : double relative_truncation_error(const DataVector& power_monitor,
92 : size_t num_modes_to_use);
93 : /// @}
94 :
95 : /*!
96 : * \brief The relative truncation error in each logical direction of the grid
97 : *
98 : * This overload is intended for visualization purposes only. It takes a tensor
99 : * component as input, so it can be used as a kernel to post-process volume data
100 : * with Python bindings (see `TransformVolumeData.py`).
101 : */
102 : template <typename VectorType, size_t Dim>
103 1 : std::array<double, Dim> relative_truncation_error(
104 : const VectorType& tensor_component, const Mesh<Dim>& mesh);
105 :
106 : /// @{
107 : /*!
108 : * \ingroup SpectralGroup
109 : * \brief Returns an estimate of the absolute truncation error in each
110 : * dimension.
111 : *
112 : * The estimate of the numerical error is given by
113 : *
114 : * \f{align*}{
115 : * \mathcal{E}\left[P_k\right] = u_\mathrm{max} \times 10^{- \mathcal{T}[P_k]},
116 : * \f}
117 : *
118 : * where \f$ u_\mathrm{max} = \mathrm{max} |u|\f$ in the corresponding element
119 : * and \f$ \mathcal{T}[P_k] \f$ is the relative error estimate
120 : * computed from the power monitors \f$ P_k \f$.
121 : *
122 : * \warning This estimate is intended for visualization purposes only.
123 : */
124 : template <typename VectorType, size_t Dim>
125 1 : std::array<double, Dim> absolute_truncation_error(
126 : const VectorType& tensor_component, const Mesh<Dim>& mesh);
127 : /// @}
128 :
129 : /// Holds convergence rate and pile up modes of a power monitor
130 1 : struct ConvergenceInfo {
131 0 : double convergence_rate{std::numeric_limits<double>::signaling_NaN()};
132 0 : double number_of_pile_up_modes{std::numeric_limits<double>::signaling_NaN()};
133 : };
134 :
135 : /*!
136 : * \ingroup SpectralGroup
137 : * \brief Returns the convergence rate and the number of pile up modes of a
138 : * power monitor as a ConvergenceInfo.
139 : *
140 : * \details Computes the convergence rate of a power monitor as a weighted
141 : * average of slopes measured using different subsets of spectral modes
142 : * in the power monitor. Equation (53) of \cite Szilagyi2014fna gives
143 : * the convergence rate $\mathcal{C}$ in terms of a power monitor $P_k$ as
144 : * \begin{equation}
145 : * \mathcal{C}(P_k) = -\frac{\sum_{k_1=0}^2\sum_{k_2=\tilde{k}_1}^{\tilde{N}-1}
146 : * \frac{\mathcal{S}(k_1,k_2)}{\epsilon + \mathcal{E}(k_1,k_2)}}{
147 : * \sum_{k_1=0}^2\sum_{k_2=\tilde{k}_1}^{\tilde{N}-1}
148 : * \frac{1}{\epsilon + \mathcal{E}(k_1,k_2)}}.
149 : * \end{equation}
150 : * Here, $\mathcal{S}(k_1,k_2)$ is the slope of a linear regression fit of
151 : * $\log_{10}(P_k)$ with $k$ satisfying $k_1\leq k \leq k_2$,
152 : * $\mathcal{E}(k_1,k_2)$ is the error of the slope in that fit,
153 : * $\epsilon=\max\left(10^{-3}\max\left(\mathcal{E}(k_1,k_2)\right),
154 : * 10^{-15}\right)$ is a small number to avoid dividing by zero in the event the
155 : * fit errors vanish,
156 : * $\max\left(\mathcal{E}(k_1,k_2)\right)$ is the maximum fit error of each
157 : * fit whose slope is included in the summation,
158 : * $\tilde{k}_1 = \min\left(k_1+4,\tilde{N}-1\right)$,
159 : * $\tilde{N} = N-N_f$, $N$ is the number of modes in the
160 : * power monitor, and the highest $N_f$ modes are filtered. Note that the
161 : * way $\epsilon$ is defined is so that it matches SpEC's definition, while
162 : * also ensuring that it is nonzero even if the error in the slope fit is
163 : * exactly zero.
164 : *
165 : * Also computes the number of pile up modes in a power monitor. Pile up
166 : * modes are modes where the power is no longer converging at the overall
167 : * convergence rate. Following Eq. (56) of \cite Szilagyi2014fna, the number of
168 : * pile up modes $\mathcal{P}$ is defined as
169 : * \begin{equation}
170 : * \mathcal{P}(P_k) = \sum_{j=2}^{\tilde{N}-2}
171 : * \exp\left[-32\left(\frac{\tilde{\mathcal{C}}_j}
172 : * {\mathcal{C}(P_k)}\right)^2\right],
173 : * \end{equation}
174 : * where $\mathcal{C}(P_k)$ is the convergence rate of the power monitor $P_k$,
175 : * the local convergence rate $\tilde{\mathcal{C}}_j$ of mode $j$ is
176 : * \begin{equation}
177 : * \tilde{\mathcal{C}}_j = -\mathcal{S}(j,\min(\tilde{N}-1,j+4)),
178 : * \end{equation}
179 : * $\mathcal{S}(k_1,k_2)$ is the slope of a linear regression fit of
180 : * $\log_{10}(P_k)$ with $k$ satisfying $k_1\leq k \leq k_2$,
181 : * $\tilde{N} = N-N_f$, $N$ is the number of modes in the
182 : * power monitor, and the highest $N_f$ modes are filtered.
183 : * The motivation of this definition is the following: if the local
184 : * convergence rate $\tilde{\mathcal{C}}_j$ is comparable to the overall
185 : * convergence rate $\mathcal{C}(P_k)$, then the $j^{\rm th}$ term in the
186 : * summation becomes $\approx \exp(-32) \approx 10^{-14}$, while if
187 : * $\tilde{\mathcal{C}}_j \ll \mathcal{C}(P_k)$, then the $j^{\rm th}$ term in
188 : * the summation is $\approx \exp(0) = 1$. Note that the coefficient value 32
189 : * is chosen to agree with SpEC.
190 : * \note The summation goes up to $\tilde{N}-2$ so that there is it least one
191 : * larger unfiltered mode for use in computing the slope. The highest mode
192 : * used when computing the slope is the highest unfiltered mode, $\tilde{N}-1$.
193 : * Mode numbers in the power monitor are zero based. These choices are off by
194 : * one vs. Eqs. (55) and (56) of \cite Szilagyi2014fna, because those formulas
195 : * apparently assume one-based indexing.
196 : * \param power_monitor The power monitor.
197 : * \param number_of_filtered_modes How many of the highest modes of the
198 : * power monitor are filtered (default 0).
199 : */
200 1 : ConvergenceInfo convergence_rate_and_number_of_pile_up_modes(
201 : const DataVector& power_monitor, size_t number_of_filtered_modes = 0);
202 :
203 : /// @{
204 : /*!
205 : * \brief Return the radial power monitor for a tensor component on a
206 : * spherical shell.
207 : *
208 : * The mesh dimensions are assumed to be ordered `(radial, theta, phi)`. The
209 : * radial grid points are contiguous, so each angular point supplies one
210 : * radial slice to the one-dimensional modal transform.
211 : */
212 1 : void spherical_shell_radial_power_monitor(gsl::not_null<DataVector*> result,
213 : const DataVector& tensor_component,
214 : const Mesh<3>& mesh);
215 :
216 1 : DataVector spherical_shell_radial_power_monitor(
217 : const DataVector& tensor_component, const Mesh<3>& mesh);
218 : /// @}
219 :
220 : /// @{
221 : /*!
222 : * \brief Return the angular power monitor for one TensorYlm component on a
223 : * spherical shell.
224 : *
225 : * The mesh dimensions are assumed to be ordered `(radial, theta, phi)`, with
226 : * `l_max == m_max`. TensorYlm coefficients use the radial dimension as the
227 : * fastest-moving extent. As reviewed in Sec. II of \cite Boyle2023,
228 : * spin-weighted spherical harmonics with `l < |spin_weight|` vanish, so these
229 : * modes are omitted from both the sum and its normalization. Set
230 : * `zero_m_is_real` for real scalar coefficients, which have no imaginary
231 : * `m=0` coefficients in Spherepack storage.
232 : */
233 1 : void spherical_shell_angular_power_monitor(
234 : gsl::not_null<DataVector*> result, const DataVector& tensor_ylm_component,
235 : const Mesh<3>& mesh, int spin_weight, bool zero_m_is_real);
236 :
237 1 : DataVector spherical_shell_angular_power_monitor(
238 : const DataVector& tensor_ylm_component, const Mesh<3>& mesh,
239 : int spin_weight, bool zero_m_is_real);
240 : /// @}
241 : } // namespace PowerMonitors
|
