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