Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include "NumericalAlgorithms/SphericalHarmonics/TensorYlm.hpp"
7 :
8 : #include <cstddef>
9 : #include <optional>
10 :
11 : #include "DataStructures/DataVector.hpp"
12 : #include "Utilities/Gsl.hpp"
13 :
14 : namespace ylm::TensorYlm {
15 :
16 : /*!
17 : * \brief Fills a sparse matrix that does a TensorYlm filter operation.
18 : *
19 : * If the CoefficientNormalization parameter is set to Standard, then
20 : * assumes that $T^{\tilde A}_{\ell' m'}$ is stored in a
21 : * Tensor<DataVector>. Multiplying (one plus) the resulting
22 : * sparse matrix by the Tensor<DataVector> is equivalent to
23 : * evaluating the right-hand side of Eq. $(\ref{eq:Filter})$.
24 : * If the CoefficientNormalization parameter is set to Spherepack, then
25 : * assumes the matrix will multiply a Tensor<DataVector> containing
26 : * ${\breve T}^{\tilde A}_{\ell' m'}$ and the filter operation is equivalent to
27 : * Eq. $(\ref{eq:FilterSpherepack})$.
28 : *
29 : * We assume that the independent components of the Tensor<DataVector>
30 : * are stored contiguously in memory, in order of the storage_index of
31 : * the Tensor. However, we do allow for a stride. This means that we
32 : * can point to memory starting with the first element of
33 : * $T^{\tilde A}_{\ell' m'}$ and multiply that by the sparse matrix we compute
34 : * here, and get the filtered result.
35 : *
36 : * The memory layout here is different than in SpEC. In SpEC, each
37 : * tensor component is stored in separately-allocated memory, so the
38 : * SpEC equivalent of the fill_filter function fills $N^2$ sparse
39 : * matrices, where $N$ is the number of independent components of the
40 : * Tensor. The advantage of the SpEC method is that each sparse matrix
41 : * is smaller, so sorting elements into the correct order while
42 : * constructing each sparse matrix is faster (sorting is > linear in
43 : * the number of matrix elements). The disadvantage of the SpEC
44 : * method is that evaluating the filter for a single Tensor involves
45 : * $N^2$ matrix-vector multiplications, whereas here evaluating the
46 : * filter involves only one matrix-vector multiplication, which should
47 : * have more efficient memory access. It is not clear which method is
48 : * faster overall without more profiling.
49 : *
50 : * ## Explicit formulas
51 : *
52 : * The following formulas come from Klinger and Scheel, in prep.
53 : *
54 : * For rank-0 tensors, the expression is simple because there
55 : * is no change of basis, only a filter based on $\ell$.
56 : * \begin{align}
57 : * F_{l m \tilde{D}}^{\ell'' m''\tilde{A}} &=
58 : * \delta(\tilde{D},\tilde{A})\delta_{\ell \ell''}\delta_{m m''}
59 : * \left[1-
60 : * \delta(\ell_{\mathrm{cut}}^-\leq \ell \leq \ell_{\mathrm{max}}) g(\ell)
61 : * \right].
62 : * \end{align}
63 : *
64 : * For rank-1 tensors, the expression for
65 : * $F_{l m \tilde{D}}^{\ell'' m''\tilde{A}}$ is
66 : * \begin{align}
67 : * F_{l m \tilde{D}}^{\ell'' m''\tilde{A}} &=
68 : * \delta(\tilde{D},\tilde{A})\delta_{\ell \ell''}\delta_{m m''}
69 : * \nonumber \\
70 : * &- (-1)^{m+m''} (-1)^{\delta(\tilde{D},\mathbf{e}_y)}
71 : * \delta(\ell,\ell'')
72 : * \sum_{\ell'=\ell_{\mathrm{cut}}^-}^{\ell_{\mathrm{max}}+1}
73 : * \frac{2\ell'+1}{2} g(\ell') \nonumber \\
74 : * &\times \sum_{j,p,m'} k_j(\tilde{D})k_p(\tilde{A})
75 : * \left(\begin{array}{rrr}
76 : * \ell&\ell'&1\cr
77 : * -m''&m'&m_j(\tilde{D})
78 : * \end{array}\right)
79 : * \left(\begin{array}{rrr}
80 : * \ell&\ell'&1\cr
81 : * -m&m'&m_p(\tilde{A})
82 : * \end{array}\right),
83 : * \end{align}
84 : * where the 6-element "matrices"
85 : * in parentheses are Wigner 3-J symbols, and where
86 : * \begin{align}
87 : * g(\ell') &=
88 : * \left\{\begin{array}{lr}
89 : * 1-f(\ell') & \ell' \leq \ell_{\mathrm{cut}}^+,\\
90 : * 1 & \ell' > \ell_{\mathrm{cut}}^+.
91 : * \end{array}\right.
92 : * \end{align}
93 : *
94 : * For second-rank tensors, the expression for
95 : * $F_{l m \tilde{D}}^{\ell'' m''\tilde{A}}$ is
96 : * \begin{align}
97 : * F_{l m \tilde{D}}^{\ell'' m''\tilde{A}}
98 : * &=
99 : * \delta(\tilde{D},\tilde{A})\delta_{\ell \ell''}\delta_{m m''}
100 : * \nonumber \\
101 : * &-
102 : * \frac{1}{4}
103 : * (-1)^{\delta(\tilde{D}_1,\mathbf{e}_y)}
104 : * (-1)^{\delta(\tilde{D}_2,\mathbf{e}_y)}
105 : * \delta_{\ell \ell''}
106 : * \nonumber \\
107 : * &\times
108 : * \sum_{\ell'=\ell_{\mathrm{cut}}^-}^{\ell_{\mathrm{max}}+2}
109 : * (2\ell'+1) g(\ell')
110 : * \sum_{u,v,p,q,\bar{\ell},\tilde{m},\bar{m},m'}
111 : * (2 \bar{\ell}+1) S(\bar{\ell},\tilde{D})
112 : * k_u(\tilde{D}_1) k_v(\tilde{D}_2)
113 : * k_p(\tilde{A}_1) k_q(\tilde{A}_2)
114 : * \nonumber \\
115 : * &\times
116 : * \left(\begin{array}{rrr}
117 : * \ell&\ell'&\bar{\ell}\cr
118 : * -m&m'&\bar{m}
119 : * \end{array}\right)
120 : * \left(\begin{array}{rrr}
121 : * 1&1&\bar{\ell}\cr
122 : * m_p(\tilde{A}_1)&m_q(\tilde{A}_2)&-\bar{m}
123 : * \end{array}\right)
124 : * \nonumber \\
125 : * &\times
126 : * \left(\begin{array}{rrr}
127 : * \ell'&\ell&\bar{\ell}\cr
128 : * -m'&m''&\tilde{m}
129 : * \end{array}\right)
130 : * \left(\begin{array}{rrr}
131 : * 1&1&\bar{\ell}\cr
132 : * m_u(\tilde{D}_1)&m_v(\tilde{D}_2)&\tilde{m}
133 : * \end{array}\right),
134 : * \label{eq:RankTwoTransformWithCut}
135 : * \end{align}
136 : * where the factor $S(\bar{\ell},\tilde{D})$ is a symmetry factor.
137 : * For a 2nd-rank tensor with no symmetries, $S(\bar{\ell},\tilde{D})$ is
138 : * unity, but for a tensor symmetric in $(\tilde{D_1},\tilde{D_2})$
139 : * the matrix elements $F_{l m \tilde{D}}^{\ell'' m''\tilde{A}}$ are
140 : * multiplied only by tensor components with $\tilde{D_1}\geq\tilde{D_2}$
141 : * and the symmetry is accounted for by setting
142 : * \begin{align}
143 : * S(\bar{\ell},\tilde{D}) &=
144 : * (1+(-1)^{\bar{\ell}})\frac{2-\delta(\tilde{D_1},\tilde{D_2})}{2}.
145 : * \end{align}
146 : *
147 : * For third-rank tensors, the expression for
148 : * $F_{l m \tilde{D}}^{\ell'' m''\tilde{A}}$ is
149 : * \begin{align}
150 : * F_{l m \tilde{D}}^{\ell'' m''\tilde{A}}
151 : * &=
152 : * \delta(\tilde{D},\tilde{A})\delta_{\ell \ell''}\delta_{m m''}
153 : * \nonumber \\
154 : * &-
155 : * \frac{1}{8}
156 : * (-1)^{\delta(\tilde{D}_1,\mathbf{e}_y)}
157 : * (-1)^{\delta(\tilde{D}_2,\mathbf{e}_y)}
158 : * (-1)^{\delta(\tilde{D}_3,\mathbf{e}_y)}
159 : * \delta_{\ell \ell''}
160 : * \nonumber \\
161 : * &\times
162 : * \sum_{\ell'=\ell_{\mathrm{cut}}^-}^{\ell_{\mathrm{max}}+3}
163 : * (2\ell'+1) g(\ell')
164 : * \sum_{u,v,w,p,q,r,\tilde{m},\bar{\ell},\bar{m},
165 : * \hat{\ell},\hat{m},\check{m},m'}
166 : * (2 \bar{\ell}+1) S(\bar{\ell},\tilde{D})
167 : * (2 \hat{\ell}+1)
168 : * \nonumber \\
169 : * &\qquad\times
170 : * k_u(\tilde{D}_2) k_v(\tilde{D}_3) k_w(\tilde{D}_1)
171 : * k_p(\tilde{A}_2) k_q(\tilde{A}_3) k_r(\tilde{A}_1)
172 : * (-1)^{\tilde{m}-\bar{m}}
173 : * \nonumber \\
174 : * &\qquad\times
175 : * \left(\begin{array}{rrr}
176 : * 1&1&\bar{\ell}\cr
177 : * m_p(\tilde{A}_2)&m_q(\tilde{A}_3)&-\bar{m}
178 : * \end{array}\right)
179 : * \left(\begin{array}{rrr}
180 : * 1&1&\bar{\ell}\cr
181 : * m_u(\tilde{D}_2)&m_v(\tilde{D}_3)&\tilde{m}
182 : * \end{array}\right)
183 : * \nonumber \\
184 : * &\qquad\times
185 : * \left(\begin{array}{rrr}
186 : * \ell&\ell'&\hat{\ell}\cr
187 : * -m&m'&\check{m}
188 : * \end{array}\right)
189 : * \left(\begin{array}{rrr}
190 : * \ell'&\ell&\hat{\ell}\cr
191 : * -m'&m''&\hat{m}
192 : * \end{array}\right)
193 : * \nonumber \\
194 : * &\qquad\times
195 : * \left(\begin{array}{rrr}
196 : * 1&\bar{\ell}&\hat{\ell}\cr
197 : * m_r(\tilde{A}_1)&\bar{m}&-\check{m}
198 : * \end{array}\right)
199 : * \left(\begin{array}{rrr}
200 : * 1&\bar{\ell}&\hat{\ell}\cr
201 : * m_w(\tilde{D}_1)&-\tilde{m}&\hat{m}
202 : * \end{array}\right).
203 : * \label{eq:RankThreeTransformWithCut}
204 : * \end{align}
205 : * As in the rank-2 case,
206 : * $S(\bar{\ell},\tilde{D})$ is a symmetry factor that is unity
207 : * for a tensor with no symmetries and is
208 : * \begin{align}
209 : * S(\bar{\ell},\tilde{D}) &=
210 : * (1+(-1)^{\bar{\ell}})\frac{2-\delta(\tilde{D_2},\tilde{D_3})}{2}
211 : * \end{align}
212 : * for a tensor symmetric on its last two indices. We don't consider
213 : * other symmetries because we don't find them in the cases we care
214 : * about.
215 : *
216 : * \tparam TensorStructure A Tensor_detail::Structure
217 : * \tparam SparseMatrixType A sparse matrix fillable by SparseMatrixFiller
218 : *
219 : * \param matrix The sparse matrix to fill
220 : * \param ell_max The maximum ylm ell value.
221 : * \param number_of_ell_modes_to_kill How many top ell modes to set to zero.
222 : * \param half_power The half power $\sigma$ for more complicated filtering.
223 : * \param coefficient_normalization Describes the normalization of coefficients.
224 : *
225 : * If half_power is std::nullopt, implements a Heaviside filter.
226 : * Otherwise, the filter is the more complicated one described in the
227 : * TensorYlm namespace documentation, with $\sigma$ equal to
228 : * half_power and $\ell_{\mathrm{cut}}^+$ equal to $\ell_{\rm max}$
229 : * minus number_of_ell_modes_to_kill.
230 : *
231 : */
232 : template <typename TensorStructure, typename SparseMatrixType>
233 1 : void fill_filter(gsl::not_null<SparseMatrixType*> matrix, size_t ell_max,
234 : size_t number_of_ell_modes_to_kill,
235 : std::optional<size_t> half_power,
236 : CoefficientNormalization coefficient_normalization);
237 :
238 : } // namespace ylm::TensorYlm
|
