Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <optional>
8 :
9 : #include "DataStructures/SimpleSparseMatrix.hpp"
10 : #include "DataStructures/Tensor/TypeAliases.hpp"
11 : #include "NumericalAlgorithms/SphericalHarmonics/TensorYlm.hpp"
12 : #include "Utilities/Gsl.hpp"
13 :
14 : /// \cond
15 : class DataVector;
16 : template <typename TagsList>
17 : class Variables;
18 : /// \endcond
19 :
20 : namespace ylm::TensorYlm {
21 :
22 : /*!
23 : * \brief Fills a sparse matrix that does a TensorYlm filter operation.
24 : *
25 : * If the CoefficientNormalization parameter is set to Standard, then
26 : * assumes that $T^{\tilde A}_{\ell' m'}$ is stored in a
27 : * Tensor<DataVector>. Multiplying (one plus) the resulting
28 : * sparse matrix by the Tensor<DataVector> is equivalent to
29 : * evaluating the right-hand side of Eq. $(\ref{eq:Filter})$.
30 : * If the CoefficientNormalization parameter is set to Spherepack, then
31 : * assumes the matrix will multiply a Tensor<DataVector> containing
32 : * ${\breve T}^{\tilde A}_{\ell' m'}$ and the filter operation is equivalent to
33 : * Eq. $(\ref{eq:FilterSpherepack})$.
34 : *
35 : * We assume that the independent components of the Tensor<DataVector>
36 : * are stored contiguously in memory, in order of the storage_index of
37 : * the Tensor. However, we do allow for a stride. This means that we
38 : * can point to memory starting with the first element of
39 : * $T^{\tilde A}_{\ell' m'}$ and multiply that by the sparse matrix we compute
40 : * here, and get the filtered result.
41 : *
42 : * The memory layout here is different than in SpEC. In SpEC, each
43 : * tensor component is stored in separately-allocated memory, so the
44 : * SpEC equivalent of the fill_filter function fills $N^2$ sparse
45 : * matrices, where $N$ is the number of independent components of the
46 : * Tensor. The advantage of the SpEC method is that each sparse matrix
47 : * is smaller, so sorting elements into the correct order while
48 : * constructing each sparse matrix is faster (sorting is > linear in
49 : * the number of matrix elements). The disadvantage of the SpEC
50 : * method is that evaluating the filter for a single Tensor involves
51 : * $N^2$ matrix-vector multiplications, whereas here evaluating the
52 : * filter involves only one matrix-vector multiplication, which should
53 : * have more efficient memory access. It is not clear which method is
54 : * faster overall without more profiling.
55 : *
56 : * ## Explicit formulas
57 : *
58 : * The following formulas come from Klinger and Scheel, in prep.
59 : *
60 : * For rank-0 tensors, the expression is simple because there
61 : * is no change of basis, only a filter based on $\ell$.
62 : * \begin{align}
63 : * F_{l m \tilde{D}}^{\ell'' m''\tilde{A}} &=
64 : * \delta(\tilde{D},\tilde{A})\delta_{\ell \ell''}\delta_{m m''}
65 : * \left[1-
66 : * \delta(\ell_{\mathrm{cut}}^-\leq \ell \leq \ell_{\mathrm{max}}) g(\ell)
67 : * \right].
68 : * \end{align}
69 : *
70 : * For rank-1 tensors, the expression for
71 : * $F_{l m \tilde{D}}^{\ell'' m''\tilde{A}}$ is
72 : * \begin{align}
73 : * F_{l m \tilde{D}}^{\ell'' m''\tilde{A}} &=
74 : * \delta(\tilde{D},\tilde{A})\delta_{\ell \ell''}\delta_{m m''}
75 : * \nonumber \\
76 : * &- (-1)^{m+m''} (-1)^{\delta(\tilde{D},\mathbf{e}_y)}
77 : * \delta(\ell,\ell'')
78 : * \sum_{\ell'=\ell_{\mathrm{cut}}^-}^{\ell_{\mathrm{max}}+1}
79 : * \frac{2\ell'+1}{2} g(\ell') \nonumber \\
80 : * &\times \sum_{j,p,m'} k_j(\tilde{D})k_p(\tilde{A})
81 : * \left(\begin{array}{rrr}
82 : * \ell&\ell'&1\cr
83 : * -m''&m'&m_j(\tilde{D})
84 : * \end{array}\right)
85 : * \left(\begin{array}{rrr}
86 : * \ell&\ell'&1\cr
87 : * -m&m'&m_p(\tilde{A})
88 : * \end{array}\right),
89 : * \end{align}
90 : * where the 6-element "matrices"
91 : * in parentheses are Wigner 3-J symbols, and where
92 : * \begin{align}
93 : * g(\ell') &=
94 : * \left\{\begin{array}{lr}
95 : * 1-f(\ell') & \ell' \leq \ell_{\mathrm{cut}}^+,\\
96 : * 1 & \ell' > \ell_{\mathrm{cut}}^+.
97 : * \end{array}\right.
98 : * \end{align}
99 : *
100 : * For second-rank tensors, the expression for
101 : * $F_{l m \tilde{D}}^{\ell'' m''\tilde{A}}$ is
102 : * \begin{align}
103 : * F_{l m \tilde{D}}^{\ell'' m''\tilde{A}}
104 : * &=
105 : * \delta(\tilde{D},\tilde{A})\delta_{\ell \ell''}\delta_{m m''}
106 : * \nonumber \\
107 : * &-
108 : * \frac{1}{4}
109 : * (-1)^{\delta(\tilde{D}_1,\mathbf{e}_y)}
110 : * (-1)^{\delta(\tilde{D}_2,\mathbf{e}_y)}
111 : * \delta_{\ell \ell''}
112 : * \nonumber \\
113 : * &\times
114 : * \sum_{\ell'=\ell_{\mathrm{cut}}^-}^{\ell_{\mathrm{max}}+2}
115 : * (2\ell'+1) g(\ell')
116 : * \sum_{u,v,p,q,\bar{\ell},\tilde{m},\bar{m},m'}
117 : * (2 \bar{\ell}+1) S(\bar{\ell},\tilde{D})
118 : * k_u(\tilde{D}_1) k_v(\tilde{D}_2)
119 : * k_p(\tilde{A}_1) k_q(\tilde{A}_2)
120 : * \nonumber \\
121 : * &\times
122 : * \left(\begin{array}{rrr}
123 : * \ell&\ell'&\bar{\ell}\cr
124 : * -m&m'&\bar{m}
125 : * \end{array}\right)
126 : * \left(\begin{array}{rrr}
127 : * 1&1&\bar{\ell}\cr
128 : * m_p(\tilde{A}_1)&m_q(\tilde{A}_2)&-\bar{m}
129 : * \end{array}\right)
130 : * \nonumber \\
131 : * &\times
132 : * \left(\begin{array}{rrr}
133 : * \ell'&\ell&\bar{\ell}\cr
134 : * -m'&m''&\tilde{m}
135 : * \end{array}\right)
136 : * \left(\begin{array}{rrr}
137 : * 1&1&\bar{\ell}\cr
138 : * m_u(\tilde{D}_1)&m_v(\tilde{D}_2)&\tilde{m}
139 : * \end{array}\right),
140 : * \label{eq:RankTwoTransformWithCut}
141 : * \end{align}
142 : * where the factor $S(\bar{\ell},\tilde{D})$ is a symmetry factor.
143 : * For a 2nd-rank tensor with no symmetries, $S(\bar{\ell},\tilde{D})$ is
144 : * unity, but for a tensor symmetric in $(\tilde{D_1},\tilde{D_2})$
145 : * the matrix elements $F_{l m \tilde{D}}^{\ell'' m''\tilde{A}}$ are
146 : * multiplied only by tensor components with $\tilde{D_1}\geq\tilde{D_2}$
147 : * and the symmetry is accounted for by setting
148 : * \begin{align}
149 : * S(\bar{\ell},\tilde{D}) &=
150 : * (1+(-1)^{\bar{\ell}})\frac{2-\delta(\tilde{D_1},\tilde{D_2})}{2}.
151 : * \end{align}
152 : *
153 : * For third-rank tensors, the expression for
154 : * $F_{l m \tilde{D}}^{\ell'' m''\tilde{A}}$ is
155 : * \begin{align}
156 : * F_{l m \tilde{D}}^{\ell'' m''\tilde{A}}
157 : * &=
158 : * \delta(\tilde{D},\tilde{A})\delta_{\ell \ell''}\delta_{m m''}
159 : * \nonumber \\
160 : * &-
161 : * \frac{1}{8}
162 : * (-1)^{\delta(\tilde{D}_1,\mathbf{e}_y)}
163 : * (-1)^{\delta(\tilde{D}_2,\mathbf{e}_y)}
164 : * (-1)^{\delta(\tilde{D}_3,\mathbf{e}_y)}
165 : * \delta_{\ell \ell''}
166 : * \nonumber \\
167 : * &\times
168 : * \sum_{\ell'=\ell_{\mathrm{cut}}^-}^{\ell_{\mathrm{max}}+3}
169 : * (2\ell'+1) g(\ell')
170 : * \sum_{u,v,w,p,q,r,\tilde{m},\bar{\ell},\bar{m},
171 : * \hat{\ell},\hat{m},\check{m},m'}
172 : * (2 \bar{\ell}+1) S(\bar{\ell},\tilde{D})
173 : * (2 \hat{\ell}+1)
174 : * \nonumber \\
175 : * &\qquad\times
176 : * k_u(\tilde{D}_2) k_v(\tilde{D}_3) k_w(\tilde{D}_1)
177 : * k_p(\tilde{A}_2) k_q(\tilde{A}_3) k_r(\tilde{A}_1)
178 : * (-1)^{\tilde{m}-\bar{m}}
179 : * \nonumber \\
180 : * &\qquad\times
181 : * \left(\begin{array}{rrr}
182 : * 1&1&\bar{\ell}\cr
183 : * m_p(\tilde{A}_2)&m_q(\tilde{A}_3)&-\bar{m}
184 : * \end{array}\right)
185 : * \left(\begin{array}{rrr}
186 : * 1&1&\bar{\ell}\cr
187 : * m_u(\tilde{D}_2)&m_v(\tilde{D}_3)&\tilde{m}
188 : * \end{array}\right)
189 : * \nonumber \\
190 : * &\qquad\times
191 : * \left(\begin{array}{rrr}
192 : * \ell&\ell'&\hat{\ell}\cr
193 : * -m&m'&\check{m}
194 : * \end{array}\right)
195 : * \left(\begin{array}{rrr}
196 : * \ell'&\ell&\hat{\ell}\cr
197 : * -m'&m''&\hat{m}
198 : * \end{array}\right)
199 : * \nonumber \\
200 : * &\qquad\times
201 : * \left(\begin{array}{rrr}
202 : * 1&\bar{\ell}&\hat{\ell}\cr
203 : * m_r(\tilde{A}_1)&\bar{m}&-\check{m}
204 : * \end{array}\right)
205 : * \left(\begin{array}{rrr}
206 : * 1&\bar{\ell}&\hat{\ell}\cr
207 : * m_w(\tilde{D}_1)&-\tilde{m}&\hat{m}
208 : * \end{array}\right).
209 : * \label{eq:RankThreeTransformWithCut}
210 : * \end{align}
211 : * As in the rank-2 case,
212 : * $S(\bar{\ell},\tilde{D})$ is a symmetry factor that is unity
213 : * for a tensor with no symmetries and is
214 : * \begin{align}
215 : * S(\bar{\ell},\tilde{D}) &=
216 : * (1+(-1)^{\bar{\ell}})\frac{2-\delta(\tilde{D_2},\tilde{D_3})}{2}
217 : * \end{align}
218 : * for a tensor symmetric on its last two indices. We don't consider
219 : * other symmetries because we don't find them in the cases we care
220 : * about.
221 : *
222 : * \tparam TensorStructure A Tensor_detail::Structure
223 : * \tparam SparseMatrixType A sparse matrix fillable by SparseMatrixFiller
224 : *
225 : * \param matrix The sparse matrix to fill
226 : * \param ell_max The maximum ylm ell value.
227 : * \param number_of_ell_modes_to_kill How many top ell modes to set to zero.
228 : * \param half_power The half power $\sigma$ for more complicated filtering.
229 : * \param coefficient_normalization Describes the normalization of coefficients.
230 : *
231 : * If half_power is std::nullopt, implements a Heaviside filter.
232 : * Otherwise, the filter is the more complicated one described in the
233 : * TensorYlm namespace documentation, with $\sigma$ equal to
234 : * half_power and $\ell_{\mathrm{cut}}^+$ equal to $\ell_{\rm max}$
235 : * minus number_of_ell_modes_to_kill.
236 : *
237 : */
238 : template <typename TensorStructure, typename SparseMatrixType>
239 1 : void fill_filter(gsl::not_null<SparseMatrixType*> matrix, size_t ell_max,
240 : size_t number_of_ell_modes_to_kill,
241 : std::optional<size_t> half_power,
242 : CoefficientNormalization coefficient_normalization);
243 :
244 : /*!
245 : * \brief Holds the filtering matrices for different matrix ranks.
246 : *
247 : * This is used for caching the matrices that are needed for a particular
248 : * system.
249 : *
250 : * Since the data being cached depends on parameters, those parameters are also
251 : * cached.
252 : */
253 1 : struct FilterMatrixHolder {
254 0 : size_t number_of_ell_modes_to_kill{};
255 0 : CoefficientNormalization coefficient_normalization{};
256 0 : std::optional<size_t> half_power;
257 :
258 0 : std::optional<SimpleSparseMatrix> scalar;
259 0 : std::optional<SimpleSparseMatrix> i;
260 0 : std::optional<SimpleSparseMatrix> ii;
261 0 : std::optional<SimpleSparseMatrix> ij;
262 0 : std::optional<SimpleSparseMatrix> kii;
263 : };
264 :
265 : /*!
266 : * \brief Fill the FilterMatrixHolder for the specific TagList (system).
267 : */
268 : template <typename TagList>
269 1 : void fill_tensor_ylm_filters(
270 : gsl::not_null<FilterMatrixHolder*> matrix, size_t ell_max,
271 : size_t number_of_ell_modes_to_kill, std::optional<size_t> half_power,
272 : CoefficientNormalization coefficient_normalization);
273 :
274 : /*!
275 : * \brief Applies TensorYlm filter in place to variables.
276 : *
277 : * When \p radial_extents is 1, `vars` and `temp_storage` are assumed to
278 : * be defined on a spherical slice, with number of grid points
279 : * corresponding to a spherical-harmonic grid of `ell_max`, and the
280 : * filter happens only on that slice.
281 : *
282 : * When \p radial_extents is > 1, `vars` and `temp_storage` are assumed to
283 : * be defined on a spherical shell of topology I1 x S2. The filter
284 : * happens in the entire volume, internally iterating over each
285 : * spherical slice at a time.
286 : *
287 : * For performance reasons, `apply_tensor_ylm_filter` does not allocate
288 : * or deallocate memory, but takes a `temp_storage` buffer. The size of
289 : * `temp_storage` should at least `radial_extents*spectral_size*num_components`,
290 : * where `num_components` is the total number of independent components in the
291 : * variable list (e.g. 5 for curved scalar wave and 50 for generalized
292 : * harmonic), and `spectral_size` is the size of the S2 Spherepack spectral
293 : * coefficient array for `ell_max`, as obtained from the member function
294 : * ylm::Spherepack::spectral_size(). Note that for S2 on Spherepack, the number
295 : * of collocation points is different than the number of spectral coefficients,
296 : * and both are different than the size of the Spherepack storage array.
297 : *
298 : * \param vars Variables at collocation points to filter.
299 : * \param temp_storage Temporary storage for the variables,
300 : * allocated outside apply_tensor_ylm_filter. See above for size requirements.
301 : * \param jac_inertial_to_grid Jacobian taking V_x from inertial to grid.
302 : * \param jac_grid_to_inertial Jacobian taking V_x from grid to inertial.
303 : * \param filter_matrices The bundle of filter matrices computed by
304 : * fill_filter. The members used depend on the ranks present in `TagList`.
305 : * \param ell_max The maximum ylm ell.
306 : * \param radial_extents The number of radial grid points, can be 1 for slices.
307 : */
308 : template <typename TagList>
309 1 : void apply_tensor_ylm_filter(
310 : gsl::not_null<Variables<TagList>*> vars,
311 : gsl::not_null<Variables<TagList>*> temp_storage,
312 : const InverseJacobian<DataVector, 3, Frame::Inertial, Frame::Grid>&
313 : jac_inertial_to_grid,
314 : const InverseJacobian<DataVector, 3, Frame::Grid, Frame::Inertial>&
315 : jac_grid_to_inertial,
316 : const FilterMatrixHolder& filter_matrices, size_t ell_max,
317 : size_t radial_extents);
318 : } // namespace ylm::TensorYlm
|
