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 :
10 : #include "Utilities/Gsl.hpp"
11 :
12 : namespace ylm::TensorYlm {
13 :
14 : /*!
15 : * \brief Fills a sparse matrix that does a TensorYlm Cartesian
16 : * to Spherical operation.
17 : *
18 : * Assumes that the input, $T^{\tilde A}_{\ell' m'}$, is stored in a
19 : * Tensor<DataVector>. Multiplying the resulting
20 : * sparse matrix by the Tensor<DataVector> is equivalent to
21 : * evaluating the right-hand side of Eq. $(\ref{eq:C2S})$.
22 : *
23 : * We assume that the independent components of the Tensor<DataVector>
24 : * are stored contiguously in memory, in order of the `storage_index` of
25 : * the Tensor. However, we do allow for a stride. This means that we
26 : * can point to memory starting with the first element of
27 : * $T^{\tilde A}_{\ell' m'}$ and multiply that by the sparse matrix we compute
28 : * here, and get the result.
29 : *
30 : * The memory layout here is different than in SpEC. In SpEC, each
31 : * tensor component is stored in separately-allocated memory, so the
32 : * SpEC equivalent of the fill_cart_to_sphere function fills $N^2$ sparse
33 : * matrices, where $N$ is the number of independent components of the
34 : * Tensor. The advantage of the SpEC method is that each sparse matrix
35 : * is smaller, so sorting elements into the correct order while
36 : * constructing each sparse matrix is faster (sorting is > linear in
37 : * the number of matrix elements). The disadvantage of the SpEC
38 : * method is that evaluating the coefficients for a single Tensor involves
39 : * $N^2$ matrix-vector multiplications, whereas here evaluating the
40 : * coefficients involves only one matrix-vector multiplication, which should
41 : * have more efficient memory access. It is not clear which method is
42 : * faster overall without more profiling.
43 : *
44 : * ## Explicit formulas
45 : *
46 : * The following formulas come from Klinger and Scheel, in prep.
47 : *
48 : * For rank-1 tensors, the expression for
49 : * $C_{\ell' m'\tilde{D}}^{\ell'' m'' B}$ is
50 : * \begin{align}
51 : * C_{\ell' m'\tilde{D}}^{\ell'' m'' B} &=
52 : * \sqrt{\frac{(2\ell'+1)(2\ell''+1)}{2}}
53 : * (-1)^{m''} (-1)^{\delta(\tilde{D},\mathbf{e}_y)} (-1)^{\delta(s_B,-1)}
54 : * \left(\begin{array}{ccc}
55 : * \ell' & \ell'' & 1 \cr
56 : * s_B & 0 & -s_B
57 : * \end{array}\right)
58 : * \nonumber \\
59 : * &\times \sum_j k_j(\tilde{D})
60 : * \left(\begin{array}{ccc}
61 : * \ell' & \ell'' & 1 \cr
62 : * -m' & m'' & -m_j(\tilde{D})
63 : * \end{array}\right),
64 : * \end{align}
65 : * where the 6-element "matrices"
66 : * in parentheses are Wigner 3-J symbols.
67 : *
68 : * For second-rank tensors, the expression for
69 : * $C_{\ell' m'\tilde{D}}^{\ell'' m'' B}$ is
70 : * \begin{align}
71 : * C_{\ell' m'\tilde{D}}^{\ell'' m'' B}
72 : * &=
73 : * \frac{1}{2}(-1)^{\delta(s_{B_1},-1)}(-1)^{\delta(s_{B_2},-1)}(-1)^{m'}
74 : * (-1)^{\delta(\tilde{D}_1,\mathbf{e}_y)}
75 : * (-1)^{\delta(\tilde{D}_2,\mathbf{e}_y)}
76 : * \sqrt{(2\ell''+1)(2\ell'+1)}
77 : * \nonumber \\
78 : * &\times
79 : * \sum_{u,v,\tilde{\ell},\tilde{m},\tilde{s}}
80 : * (2 \tilde{\ell}+1) (-1)^{\tilde{s}} S(\tilde{\ell},\tilde{D})
81 : * k_u(\tilde{D}_1) k_v(\tilde{D}_2)
82 : * \left(\begin{array}{ccc}
83 : * \ell' & \ell'' & \tilde{\ell} \cr
84 : * -m' & m'' & \tilde{m}
85 : * \end{array}\right)
86 : * \nonumber \\
87 : * &\qquad\times
88 : * \left(\begin{array}{ccc}
89 : * \ell' & \ell'' & \tilde{\ell} \cr
90 : * s_{B_1}+s_{B_2} & 0 & -\tilde{s}
91 : * \end{array}\right)
92 : * \left(\begin{array}{ccc}
93 : * 1 & 1 & \tilde{\ell} \cr
94 : * -m_u(\tilde{D}_1)&-m_v(\tilde{D}_2)&-\tilde{m}
95 : * \end{array}\right)
96 : * \left(\begin{array}{ccc}
97 : * 1 & 1 & \tilde{\ell} \cr
98 : * -s_{B_1}&-s_{B_2}&\tilde{s}
99 : * \end{array}\right),
100 : * \end{align}
101 : * where the factor $S(\tilde{\ell},\tilde{D})$ is a symmetry factor.
102 : * For a 2nd-rank tensor with no symmetries, $S(\tilde{\ell},\tilde{D})$ is
103 : * unity, but for a tensor symmetric in $(\tilde{D_1},\tilde{D_2})$
104 : * the matrix elements $C_{\ell' m'\tilde{D}}^{\ell'' m'' B}$ are
105 : * multiplied only by tensor components with $\tilde{D_1}\geq\tilde{D_2}$
106 : * and the symmetry is accounted for by setting
107 : * \begin{align}
108 : * S(\tilde{\ell},\tilde{D}) &=
109 : * (1+(-1)^{\tilde{\ell}})\frac{2-\delta(\tilde{D_1},\tilde{D_2})}{2}.
110 : * \end{align}
111 : *
112 : * For third-rank tensors, the expression for
113 : * $C_{\ell' m'\tilde{D}}^{\ell'' m'' B}$ is
114 : * \begin{align}
115 : * C_{\ell' m'\tilde{D}}^{\ell'' m'' B}
116 : * &=
117 : * (-1)^{m'}
118 : * (-1)^{\delta(s_{B_1},-1)}
119 : * (-1)^{\delta(s_{B_2},-1)}
120 : * (-1)^{\delta(s_{B_3},-1)}
121 : * (-1)^{\delta(\tilde{D}_1,\mathbf{e}_y)}
122 : * (-1)^{\delta(\tilde{D}_2,\mathbf{e}_y)}
123 : * (-1)^{\delta(\tilde{D}_3,\mathbf{e}_y)}
124 : * \nonumber \\
125 : * &\times
126 : * \sqrt{\frac{(2\ell''+1)(2\ell'+1)}{8}}
127 : * \sum_{u,v,w,\tilde{\ell},\tilde{m},\tilde{s},\hat{\ell},\hat{m},\hat{s}}
128 : * (2 \tilde{\ell}+1) (2 \hat{\ell}+1) (-1)^{\tilde{s}+\hat{s}+\tilde{m}}
129 : * k_u(\tilde{D}_2) k_v(\tilde{D}_3) k_w(\tilde{D}_1)
130 : * \nonumber \\
131 : * &\qquad\times
132 : * S(\tilde{\ell},\tilde{D})
133 : * \left(\begin{array}{ccc}
134 : * \ell' & \ell'' & \hat{\ell} \cr
135 : * -m' & m'' & \hat{m}
136 : * \end{array}\right)
137 : * \left(\begin{array}{ccc}
138 : * \ell' & \ell'' & \hat{\ell} \cr
139 : * s_{B_1}+s_{B_2}+s_{B_3} & 0 & -\hat{s}
140 : * \end{array}\right)
141 : * \nonumber \\
142 : * &\qquad\times
143 : * \left(\begin{array}{ccc}
144 : * 1 & 1 & \tilde{\ell} \cr
145 : * -m_u(\tilde{D}_2)&-m_v(\tilde{D}_3)&-\tilde{m}
146 : * \end{array}\right)
147 : * \left(\begin{array}{ccc}
148 : * 1 & 1 & \tilde{\ell} \cr
149 : * -s_{B_2}&-s_{B_3}&\tilde{s}
150 : * \end{array}\right)
151 : * \nonumber \\
152 : * &\qquad\times
153 : * \left(\begin{array}{ccc}
154 : * 1 & \tilde{\ell} & \hat{\ell} \cr
155 : * -m_w(\tilde{D}_1)&\tilde{m}&-\hat{m}
156 : * \end{array}\right)
157 : * \left(\begin{array}{ccc}
158 : * 1 & \tilde{\ell} & \hat{\ell} \cr
159 : * -s_{B_1}&-\tilde{s}&\hat{s}
160 : * \end{array}\right).
161 : * \end{align}
162 : * As in the rank-2 case,
163 : * $S(\tilde{\ell},\tilde{D})$ is a symmetry factor that is unity
164 : * for a tensor with no symmetries and is
165 : * \begin{align}
166 : * S(\tilde{\ell},\tilde{D}) &=
167 : * (1+(-1)^{\tilde{\ell}})\frac{2-\delta(\tilde{D_2},\tilde{D_3})}{2}
168 : * \end{align}
169 : * for a tensor symmetric on its last two indices. We don't consider
170 : * other symmetries because we don't find them in the cases we care
171 : * about.
172 : *
173 : * \tparam TensorStructure A Tensor_detail::Structure
174 : * \tparam SparseMatrixType A sparse matrix fillable by SparseMatrixFiller
175 : *
176 : * \param matrix The sparse matrix to fill
177 : * \param ell_max The maximum ylm ell value
178 : *
179 : */
180 : template <typename TensorStructure, typename SparseMatrixType>
181 1 : void fill_cart_to_sphere(gsl::not_null<SparseMatrixType*> matrix,
182 : size_t ell_max);
183 :
184 : } // namespace ylm::TensorYlm
|
