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