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^B_{\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:S2C})$.
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^B_{\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_sphere_to_cart 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{A}}_{\ell m B}$ is
50 : * \begin{align}
51 : * C^{\ell' m' \tilde{A}}_{\ell m B} &=
52 : * (-1)^{\delta(s_B,-1)-m}\sqrt{\frac{(2 \ell+1)(2 \ell' + 1)}{2}}
53 : * \sum_j k_j(\tilde{A})
54 : * \left(\begin{array}{ccc}
55 : * \ell & \ell' & 1 \cr
56 : * 0 & -s_B & s_B
57 : * \end{array}\right)
58 : * \left(\begin{array}{ccc}
59 : * \ell & \ell' & 1 \cr
60 : * -m & m' & m_j(\tilde{A})
61 : * \end{array}\right),
62 : * \end{align}
63 : * where the 6-element "matrices"
64 : * in parentheses are Wigner 3-J symbols.
65 : *
66 : * For second-rank tensors, the expression for
67 : * $C^{\ell' m' \tilde{A}}_{\ell m B}$ is
68 : * \begin{align}
69 : * C^{\ell' m' \tilde{A}}_{\ell m B}
70 : * &=
71 : * \frac{1}{2}(-1)^{\delta(s_{B_1},-1)}(-1)^{\delta(s_{B_2},-1)}(-1)^{m'}
72 : * \sqrt{(2\ell+1)(2\ell'+1)}
73 : * \nonumber \\
74 : * &\times
75 : * \sum_{p,q,\bar{\ell},\bar{m},\bar{s}}
76 : * (2 \bar{\ell}+1) (-1)^{\bar{s}}
77 : * k_p(\tilde{A}_1) k_q(\tilde{A}_2) S(\bar{\ell},B)
78 : * \left(\begin{array}{ccc}
79 : * \ell & \ell' & \bar{\ell} \cr
80 : * 0&-s_{B_1}-s_{B_2}&-\bar{s}
81 : * \end{array}\right)
82 : * \nonumber \\
83 : * &\qquad\times
84 : * \left(\begin{array}{ccc}
85 : * \ell & \ell' & \bar{\ell} \cr
86 : * -m & m' & \bar{m}
87 : * \end{array}\right)
88 : * \left(\begin{array}{ccc}
89 : * 1 & 1 & \bar{\ell} \cr
90 : * m_p(\tilde{A}_1)&m_q(\tilde{A}_2)&-\bar{m}
91 : * \end{array}\right)
92 : * \left(\begin{array}{ccc}
93 : * 1 & 1 & \bar{\ell} \cr
94 : * s_{B_1}&s_{B_2}&\bar{s}
95 : * \end{array}\right),
96 : * \end{align}
97 : * where the factor $S(\bar{\ell},B)$ is a symmetry factor.
98 : * For a 2nd-rank tensor with no symmetries, $S(\bar{\ell},B)$ is
99 : * unity, but for a tensor symmetric in $(B_1,B_2)$
100 : * the matrix elements $C^{\ell' m' \tilde{A}}_{\ell m B}$
101 : * are multiplied only by tensor components with $B_1\geq B_2$
102 : * and the symmetry is accounted for by setting
103 : * \begin{align}
104 : * S(\bar{\ell},B) &=
105 : * (1+(-1)^{\bar{\ell}})\frac{2-\delta(B_1,B_2)}{2}.
106 : * \end{align}
107 : *
108 : * For third-rank tensors, the expression for
109 : * $C^{\ell' m' \tilde{A}}_{\ell m B}$ is
110 : * \begin{align}
111 : * C^{\ell' m' \tilde{A}}_{\ell m B}
112 : * &=
113 : * (-1)^{m'}
114 : * (-1)^{\delta(s_{B_1},-1)}
115 : * (-1)^{\delta(s_{B_2},-1)}
116 : * (-1)^{\delta(s_{B_3},-1)}
117 : * \sqrt{\frac{(2\ell+1)(2\ell'+1)}{8}}
118 : * \nonumber \\
119 : * &\times
120 : * \sum_{p,q,r,\bar{\ell},\bar{m},\bar{s},\check{\ell},\check{m},\check{s}}
121 : * (2 \bar{\ell}+1) (2 \check{\ell}+1) (-1)^{\bar{s}+\check{s}-\bar{m}}
122 : * k_p(\tilde{A}_2) k_q(\tilde{A}_3) k_r(\tilde{A}_1) S(\bar{\ell},B)
123 : * \nonumber \\
124 : * &\qquad\times
125 : * \left(\begin{array}{ccc}
126 : * \ell & \ell' & \check{\ell} \cr
127 : * -m & m' & \check{m}
128 : * \end{array}\right)
129 : * \left(\begin{array}{ccc}
130 : * \ell & \ell' & \check{\ell} \cr
131 : * 0&-s_{B_1}-s_{B_2}-s_{B_3}&-\check{s}
132 : * \end{array}\right)
133 : * \nonumber \\
134 : * &\qquad\times
135 : * \left(\begin{array}{ccc}
136 : * 1 & 1 & \bar{\ell} \cr
137 : * m_p(\tilde{A}_2)&m_q(\tilde{A}_3)&-\bar{m}
138 : * \end{array}\right)
139 : * \left(\begin{array}{ccc}
140 : * 1 & 1 & \bar{\ell} \cr
141 : * s_{B_2}&s_{B_3}&\bar{s}
142 : * \end{array}\right)
143 : * \nonumber \\
144 : * &\qquad\times
145 : * \left(\begin{array}{ccc}
146 : * 1 & \bar{\ell} & \check{\ell} \cr
147 : * m_r(\tilde{A}_1)&\bar{m}&-\check{m}
148 : * \end{array}\right)
149 : * \left(\begin{array}{ccc}
150 : * 1 & \bar{\ell} & \check{\ell} \cr
151 : * s_{B_1}&-\bar{s}&\check{s}
152 : * \end{array}\right),
153 : * \end{align}
154 : * As in the rank-2 case,
155 : * $S(\bar{\ell},B)$ is a symmetry factor that is unity
156 : * for a tensor with no symmetries and is
157 : * \begin{align}
158 : * S(\bar{\ell},B) &=
159 : * (1+(-1)^{\bar{\ell}})\frac{2-\delta(B_2,B_3)}{2}.
160 : * \end{align}
161 : * for a tensor symmetric on its last two indices. We don't consider
162 : * other symmetries because we don't find them in the cases we care
163 : * about.
164 : *
165 : * \tparam TensorStructure A Tensor_detail::Structure
166 : * \tparam SparseMatrixType A sparse matrix fillable by SparseMatrixFiller
167 : *
168 : * \param matrix The sparse matrix to fill
169 : * \param ell_max The maximum ylm ell value
170 : *
171 : */
172 : template <typename TensorStructure, typename SparseMatrixType>
173 1 : void fill_sphere_to_cart(gsl::not_null<SparseMatrixType*> matrix,
174 : size_t ell_max);
175 :
176 : } // namespace ylm::TensorYlm
|
