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