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 : #include <vector>
10 :
11 : #include "DataStructures/DataVector.hpp"
12 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
13 : #include "NumericalAlgorithms/TensorYlm/Helpers.hpp"
14 : #include "Utilities/ConstantExpressions.hpp"
15 : #include "Utilities/Gsl.hpp"
16 : #include "Utilities/TMPL.hpp"
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 :
242 : /*!
243 : * \brief Return the RMS radial power monitor across all components of a
244 : * tensor on a spherical shell.
245 : *
246 : * Combines `spherical_shell_radial_power_monitor` for each of
247 : * `tensor.size()` components in quadrature, normalized by the number of
248 : * components.
249 : */
250 : template <typename TensorType>
251 1 : DataVector spherical_shell_tensor_radial_power_monitor(const TensorType& tensor,
252 : const Mesh<3>& mesh) {
253 : DataVector squared_power(mesh.extents(0), 0.0);
254 : DataVector component_power{};
255 : for (size_t component = 0; component < tensor.size(); ++component) {
256 : spherical_shell_radial_power_monitor(make_not_null(&component_power),
257 : tensor[component], mesh);
258 : squared_power += square(component_power);
259 : }
260 : squared_power = sqrt(squared_power / static_cast<double>(tensor.size()));
261 : return squared_power;
262 : }
263 :
264 : /*!
265 : * \brief Number of independent TensorYlm coefficients contributing to
266 : * angular degree `ell`, summed over `radial_extents` radial points.
267 : *
268 : * Returns 0 when `ell < |spin_weight|`, since spin-weighted spherical
269 : * harmonics vanish there (such terms are then excluded from
270 : * `accumulate_spherical_shell_tensor_angular_power()`'s sum and normalization).
271 : * See `SpherepackIterator` for the meaning of `zero_m_is_real`.
272 : */
273 1 : size_t spherical_shell_number_of_angular_coefficients(size_t ell,
274 : int spin_weight,
275 : bool zero_m_is_real,
276 : size_t radial_extents);
277 :
278 : /*!
279 : * \brief Accumulate weighted-squared angular TensorYlm power and mode counts
280 : * for all components of a tensor on a spherical shell.
281 : *
282 : * Adds to `weighted_squared_power` and `counts` in place, so this can be
283 : * called repeatedly to combine several tensors into the same angular power
284 : * monitor before calling `normalize_spherical_shell_angular_power`.
285 : */
286 : template <typename TensorType>
287 1 : void accumulate_spherical_shell_tensor_angular_power(
288 : const gsl::not_null<DataVector*> weighted_squared_power,
289 : const gsl::not_null<std::vector<size_t>*> counts, const TensorType& tensor,
290 : const Mesh<3>& mesh) {
291 : const size_t radial_extents = mesh.extents(0);
292 : const size_t ell_max = mesh.extents(1) - 1;
293 : constexpr bool zero_m_is_real = TensorType::rank() == 0;
294 : DataVector component_power{};
295 : for (size_t component = 0; component < tensor.size(); ++component) {
296 : const int spin_weight = ylm::TensorYlm::helpers::component_spin_weight<
297 : typename TensorType::structure>(component);
298 : spherical_shell_angular_power_monitor(make_not_null(&component_power),
299 : tensor[component], mesh, spin_weight,
300 : zero_m_is_real);
301 : for (size_t ell = 0; ell <= ell_max; ++ell) {
302 : const size_t component_count =
303 : spherical_shell_number_of_angular_coefficients(
304 : ell, spin_weight, zero_m_is_real, radial_extents);
305 : (*weighted_squared_power)[ell] +=
306 : static_cast<double>(component_count) * square(component_power[ell]);
307 : (*counts)[ell] += component_count;
308 : }
309 : }
310 : }
311 :
312 : /*!
313 : * \brief Normalize an angular power monitor accumulator in place.
314 : *
315 : * Sets `power[ell] = sqrt(power[ell] / counts[ell])`, or 0 when
316 : * `counts[ell] == 0`.
317 : */
318 1 : void normalize_spherical_shell_angular_power(gsl::not_null<DataVector*> power,
319 : const std::vector<size_t>& counts);
320 :
321 : /*!
322 : * \brief Accumulate squared ZernikeB3 Jacobi spectral coefficients for one
323 : * TensorYlm component into radial and angular power bins.
324 : *
325 : * `spec_buf` has layout `spec_buf[s * n_r + i_r]` where `s` is the SPHEREPACK
326 : * offset and `i_r` is the radial collocation index. Modes are binned radially
327 : * via `radial_mode = (n_total + 1) / 2` where `n_total = l + 2 * k_spec`, and
328 : * angularly by degree \f$\ell\f$. Modes with `l < |spin_weight|` are skipped.
329 : *
330 : * `offsets_by_l` must be pre-computed for the correct `zero_m_is_real` value
331 : * of this component. `gathered` and `modal_buf` are caller-owned scratch
332 : * buffers of size `max_n_modes_l * n_r` each, where
333 : * `max_n_modes_l = 2 * (l_max + 1)`.
334 : */
335 1 : void accumulate_b3_tensor_component_sums(
336 : gsl::not_null<DataVector*> sum_sq_radial,
337 : gsl::not_null<DataVector*> counts_radial,
338 : gsl::not_null<DataVector*> sum_sq_angular,
339 : gsl::not_null<DataVector*> counts_angular, const double* spec_buf,
340 : size_t n_r, size_t n_r_max, int spin_weight,
341 : const std::vector<std::vector<size_t>>& offsets_by_l, double* gathered,
342 : double* modal_buf);
343 :
344 : /*!
345 : * \brief Accumulate squared ZernikeB3 spectral coefficients for all components
346 : * of a TensorYlm tensor into radial and angular power bins.
347 : *
348 : * Selects `offsets_by_l_real` for rank-0 (scalar) tensors and
349 : * `offsets_by_l_complex` for all higher-rank tensors. The spin weight for each
350 : * component is determined from the tensor structure.
351 : */
352 : template <typename TensorType>
353 1 : void accumulate_b3_tensor_sums(
354 : gsl::not_null<DataVector*> sum_sq_radial,
355 : gsl::not_null<DataVector*> counts_radial,
356 : gsl::not_null<DataVector*> sum_sq_angular,
357 : gsl::not_null<DataVector*> counts_angular, const TensorType& tensor,
358 : size_t n_r, size_t n_r_max,
359 : const std::vector<std::vector<size_t>>& offsets_by_l_real,
360 : const std::vector<std::vector<size_t>>& offsets_by_l_complex,
361 : double* gathered, double* modal_buf) {
362 : constexpr bool zero_m_is_real = TensorType::rank() == 0;
363 : const auto& offsets_by_l =
364 : zero_m_is_real ? offsets_by_l_real : offsets_by_l_complex;
365 : for (size_t component = 0; component < tensor.size(); ++component) {
366 : const int spin_weight = ylm::TensorYlm::helpers::component_spin_weight<
367 : typename TensorType::structure>(component);
368 : accumulate_b3_tensor_component_sums(
369 : sum_sq_radial, counts_radial, sum_sq_angular, counts_angular,
370 : tensor[component].data(), n_r, n_r_max, spin_weight, offsets_by_l,
371 : gathered, modal_buf);
372 : }
373 : }
374 :
375 : /*!
376 : * \brief Normalize a B3 power monitor accumulator in place.
377 : *
378 : * Sets `result[i] = sqrt(sum_sq[i] / counts[i])`, or 0 when `counts[i] == 0`.
379 : */
380 1 : void normalize_b3_power(gsl::not_null<DataVector*> result,
381 : const DataVector& sum_sq, const DataVector& counts);
382 :
383 : } // namespace PowerMonitors
|
