|
SpECTRE
v2025.08.19
|
Holds functions that converts spatial tensors between scalar-Ylm and tensor-Ylm basis, and functions that filter by doing such a conversion. More...
Functions | |
| template<typename TensorStructure , typename SparseMatrixType > | |
| void | fill_cart_to_sphere (gsl::not_null< SparseMatrixType * > matrix, size_t ell_max) |
| Fills a sparse matrix that does a TensorYlm Cartesian to Spherical operation. More... | |
| template<typename TensorStructure , typename SparseMatrixType > | |
| void | fill_filter (gsl::not_null< SparseMatrixType * > matrix, size_t ell_max, size_t number_of_ell_modes_to_kill, std::optional< size_t > half_power) |
| Fills a sparse matrix that does a TensorYlm filter operation. More... | |
| template<typename TensorStructure , typename SparseMatrixType > | |
| void | fill_sphere_to_cart (gsl::not_null< SparseMatrixType * > matrix, size_t ell_max) |
| Fills a sparse matrix that does a TensorYlm Cartesian to Spherical operation. More... | |
Holds functions that converts spatial tensors between scalar-Ylm and tensor-Ylm basis, and functions that filter by doing such a conversion.
We expand a spatial tensor of arbitrary rank in terms of Cartesian components as follows:
\begin{align} {\mathbf T} &= \sum_{\ell,m,\tilde{A}} {}_0 Y_{\ell m} T^{\tilde A}_{\ell m}\mathbf{e}_{\tilde A}, \end{align}
where \({}_0 Y_{\ell m}\) are the usual scalar spherical harmonics, \(T^{\tilde A}_{\ell m}\) are complex-valued expansion coefficients, \(\tilde A\) is a multi-index \(\tilde{A}=(\tilde{a}_1,\tilde{a}_2,\tilde{a}_3,\ldots)\), where each of the \(\tilde{a}_i\) take on values from 0 to 2, referring to either \(\mathbf{e}_x\), \(\mathbf{e}_y\), or \(\mathbf{e}_z\). The sum over \(\tilde A\) goes over all \(\tilde A\) of the same rank.
Similarly, we can expand the same spatial tensor in terms of a complex tetrad:
\begin{align} {\mathbf T} &= \sum_{\ell,m,A} {}_{s(\!A\!)}Y_{\ell m} T^A_{\ell m} \mathbf{e}_A, \end{align}
where \({}_sY_{\ell m}\) are the spin-weighted harmonics, \(A\) is a multi-index \(A=(a_1,a_2,a_3,\ldots)\), where each of the \(a_i\) refer to either \(\mathbf{l}\), \(\mathbf{m}\), or \(\mathbf{\bar{m}}\), and \(s(\!A\!)\) is the spin weight, which is a function of \(A\): each \(\mathbf{l}\) in \(A\) adds the value zero, each \(\mathbf{m}\) adds the value -1, and each \(\mathbf{\bar{m}}\) adds the value +1. (These values are the spin weights of the complex conjugates of the basis vectors.)
To allow for a compact notation where we don't need to write different equations for different combinations of basis vectors, we define two auxiliary quantities \(k_j\) and \(m_j\) that are functions of the Cartesian basis vectors:
\begin{align} k_j(\mathbf{e}_x) &= -j,\\ k_j(\mathbf{e}_y) &= i,\\ k_j(\mathbf{e}_z) &= 1/\sqrt{2},\\ m_j(\mathbf{e}_x) &= j,\\ m_j(\mathbf{e}_y) &= j,\\ m_j(\mathbf{e}_z) &= 0, \end{align}
where the \(i\) above is the imaginary unit \(i=\sqrt{-1}\). These quantities will be used in explicit formulas for transformations between expansion coefficients.
We can define the following transformations between the expansion coefficients \(T^A_{\ell m}\) and \(T^{\tilde A}_{\ell m}\):
\begin{align} T^{\tilde B}_{\ell' m'} &= \sum_{\ell,m,A} C_{\ell' m' A}^{\ell m\tilde{B}} T^A_{\ell m},\\ T^{B}_{\ell' m'} &= \sum_{\ell m \tilde{A}} C_{\ell' m'\tilde{A}}^{\ell mB} T^{\tilde A}_{\ell m}, \end{align}
where the above transformations define the coefficients \(C_{\ell' m' A}^{\ell m\tilde{B}}\) and \(C_{\ell' m'\tilde{A}}^{\ell mB}\). Analytic expressions for these coefficients are derived in Klinger and Scheel (in prep) and will be used here to compute the coefficients.
For real-valued tensors, the coefficients \(T^A_{\ell m}\) and \(T^{\tilde A}_{\ell m}\) for negative \(m\) can be computed from the complex conjugates of the same coefficients with positive \(m\). Therefore we store and loop over only nonnegative \(m\). This means we can write
\begin{align} T^{\tilde B}_{\ell m} &= \sum_{A, \ell',m'\geq 0} \left[ C_{\ell m A}^{\ell' m'\tilde{B}} T^A_{\ell' m'} +{\hat C}_{\ell m A}^{\ell' m'\tilde{B}} T^A_{\ell' m'}{}^\star \right]\label{eq:S2C},\\ T^{B}_{\ell m} &= \sum_{\tilde{A},\ell',m'\geq 0} \left[ C_{\ell m\tilde{A}}^{\ell' m' B} T^{\tilde A}_{\ell' m'} +{\hat C}_{\ell m\tilde{A}}^{\ell' m' B} T^{\tilde A}_{\ell' m'}{}^\star \right]\label{eq:C2S}, \end{align}
where
\begin{align} {\hat C}_{\ell m A}^{\ell' m'\tilde{B}} &= (1-\delta_{0m'})C_{\ell m A^\star}^{\ell' -m'\tilde{B}} (-1)^{S(A)+m'},\\ {\hat C}_{\ell m\tilde{A}}^{\ell' m' B} &= (1-\delta_{0m'})C_{\ell m \tilde{A}}^{\ell' -m'B} (-1)^{m'}. \end{align}
The functions FillCartToSphere and FillSphereToCart fill sparse matrices that encode Eqs. \((\ref{eq:C2S})\) and \((\ref{eq:S2C})\).
Consider starting with coefficients \(T^{\tilde A}_{\ell' m'}\), transforming to spin-weighted harmonics, applying a filter to the spin-weighted harmonic coefficients, and transforming back. We can write this entire operation as
\begin{align} F_{\ell m \tilde{D}}^{\ell'' m''\tilde{A}} &= \sum_{\ell'=0}^{\ell_{\mathrm{cut}}^--1} C^{\ell' m' \tilde{A}}_{\ell m B} C^{\ell'' m'' B}_{\ell' m' \tilde{D}} + \sum_{\ell'=\ell_{\mathrm{cut}}^-}^{\ell_{\mathrm{cut}}^+} C^{\ell' m' \tilde{A}}_{\ell m B} C^{\ell'' m'' B}_{\ell' m' \tilde{D}} f(\ell'), \end{align}
where \(f(\ell')\) is a filter function, \(\ell_{\mathrm{cut}}^{-}\) is the smallest \(\ell\) mode that is filtered, and \(\ell_{\mathrm{cut}}^{+}\) is the largest \(\ell\) mode that is retained.
There are two cases we care about. The first is simple "Heaviside filtering", where \(\ell_{\mathrm{cut}}^{-} = \ell_{\mathrm{cut}}^{+}\) and \(f(\ell')=1\).
The second case is
\begin{align} f(\ell') &= \exp\left[-\alpha \left(\frac{\ell'}{\ell_{\mathrm{cut}}^++1}\right)^{2 \sigma}\right],\\ \ell_{\mathrm{cut}}^- &= \mathrm{ceil}\left[(\ell_{\mathrm{cut}}^++1) \left(\frac{\epsilon}{\alpha}\right)^{1/(2\sigma)}\right], \end{align}
where \(\alpha\) is a parameter we choose to be 36, \(\sigma\) is an integer parameter typically between 28 and 32, and \(\epsilon\) is machine roundoff. Note that the filter remains smooth at \(\ell' = \ell_{\mathrm{cut}}^-\) and it reduces to the simple filter as \(\sigma \to \infty\).
As with the transformations, we sum over only nonnegative \(m\), so we can write the action of the filter as
\begin{align} T^{\tilde B}_{\ell m} &= \sum_{{\tilde A}, \ell',m'\geq 0} \left[ F_{\ell m {\tilde A}}^{\ell' m'\tilde{B}} T^{\tilde A}_{\ell' m'} +{\hat F}_{\ell m {\tilde A}}^{\ell' m'\tilde{B}} T^{\tilde A}_{\ell' m'}{}^\star \right]\label{eq:Filter}, \end{align}
where
\begin{align} {\hat F}_{\ell m\tilde{A}}^{\ell' m' \tilde{B}} &= (1-\delta_{0m'})F_{\ell m \tilde{A}}^{\ell' -m' \tilde{B}} (-1)^{m'}. \end{align}
The function FillFilter encapsulates Eq. \((\ref{eq:Filter})\).
| void ylm::TensorYlm::fill_cart_to_sphere | ( | gsl::not_null< SparseMatrixType * > | matrix, |
| size_t | ell_max | ||
| ) |
Fills a sparse matrix that does a TensorYlm Cartesian to Spherical operation.
Assumes that the input, \(T^{\tilde A}_{\ell' m'}\), is stored in a Tensor<DataVector>. Multiplying the resulting sparse matrix by the Tensor<DataVector> is equivalent to evaluating the right-hand side of Eq. \((\ref{eq:C2S})\).
We assume that the independent components of the Tensor<DataVector> are stored contiguously in memory, in order of the storage_index of the Tensor. However, we do allow for a stride. This means that we can point to memory starting with the first element of \(T^{\tilde A}_{\ell' m'}\) and multiply that by the sparse matrix we compute here, and get the result.
The memory layout here is different than in SpEC. In SpEC, each tensor component is stored in separately-allocated memory, so the SpEC equivalent of the fill_cart_to_sphere function fills \(N^2\) sparse matrices, where \(N\) is the number of independent components of the Tensor. The advantage of the SpEC method is that each sparse matrix is smaller, so sorting elements into the correct order while constructing each sparse matrix is faster (sorting is > linear in the number of matrix elements). The disadvantage of the SpEC method is that evaluating the coefficients for a single Tensor involves \(N^2\) matrix-vector multiplications, whereas here evaluating the coefficients involves only one matrix-vector multiplication, which should have more efficient memory access. It is not clear which method is faster overall without more profiling.
The following formulas come from Klinger and Scheel, in prep.
For rank-1 tensors, the expression for \(C_{\ell' m'\tilde{D}}^{\ell'' m'' B}\) is
\begin{align} C_{\ell' m'\tilde{D}}^{\ell'' m'' B} &= \sqrt{\frac{(2\ell'+1)(2\ell''+1)}{2}} (-1)^{m''} (-1)^{\delta(\tilde{D},\mathbf{e}_y)} (-1)^{\delta(s_B,-1)} \left(\begin{array}{ccc} \ell' & \ell'' & 1 \cr s_B & 0 & -s_B \end{array}\right) \nonumber \\ &\times \sum_j k_j(\tilde{D}) \left(\begin{array}{ccc} \ell' & \ell'' & 1 \cr -m' & m'' & -m_j(\tilde{D}) \end{array}\right), \end{align}
where the 6-element "matrices" in parentheses are Wigner 3-J symbols.
For second-rank tensors, the expression for \(C_{\ell' m'\tilde{D}}^{\ell'' m'' B}\) is
\begin{align} C_{\ell' m'\tilde{D}}^{\ell'' m'' B} &= \frac{1}{2}(-1)^{\delta(s_{B_1},-1)}(-1)^{\delta(s_{B_2},-1)}(-1)^{m'} (-1)^{\delta(\tilde{D}_1,\mathbf{e}_y)} (-1)^{\delta(\tilde{D}_2,\mathbf{e}_y)} \sqrt{(2\ell''+1)(2\ell'+1)} \nonumber \\ &\times \sum_{u,v,\tilde{\ell},\tilde{m},\tilde{s}} (2 \tilde{\ell}+1) (-1)^{\tilde{s}} S(\tilde{\ell},\tilde{D}) k_u(\tilde{D}_1) k_v(\tilde{D}_2) \left(\begin{array}{ccc} \ell' & \ell'' & \tilde{\ell} \cr -m' & m'' & \tilde{m} \end{array}\right) \nonumber \\ &\qquad\times \left(\begin{array}{ccc} \ell' & \ell'' & \tilde{\ell} \cr s_{B_1}+s_{B_2} & 0 & -\tilde{s} \end{array}\right) \left(\begin{array}{ccc} 1 & 1 & \tilde{\ell} \cr -m_u(\tilde{D}_1)&-m_v(\tilde{D}_2)&-\tilde{m} \end{array}\right) \left(\begin{array}{ccc} 1 & 1 & \tilde{\ell} \cr -s_{B_1}&-s_{B_2}&\tilde{s} \end{array}\right), \end{align}
where the factor \(S(\tilde{\ell},\tilde{D})\) is a symmetry factor. For a 2nd-rank tensor with no symmetries, \(S(\tilde{\ell},\tilde{D})\) is unity, but for a tensor symmetric in \((\tilde{D_1},\tilde{D_2})\) the matrix elements \(C_{\ell' m'\tilde{D}}^{\ell'' m'' B}\) are multiplied only by tensor components with \(\tilde{D_1}\geq\tilde{D_2}\) and the symmetry is accounted for by setting
\begin{align} S(\tilde{\ell},\tilde{D}) &= (1+(-1)^{\tilde{\ell}})\frac{2-\delta(\tilde{D_1},\tilde{D_2})}{2}. \end{align}
For third-rank tensors, the expression for \(C_{\ell' m'\tilde{D}}^{\ell'' m'' B}\) is
\begin{align} C_{\ell' m'\tilde{D}}^{\ell'' m'' B} &= (-1)^{m'} (-1)^{\delta(s_{B_1},-1)} (-1)^{\delta(s_{B_2},-1)} (-1)^{\delta(s_{B_3},-1)} (-1)^{\delta(\tilde{D}_1,\mathbf{e}_y)} (-1)^{\delta(\tilde{D}_2,\mathbf{e}_y)} (-1)^{\delta(\tilde{D}_3,\mathbf{e}_y)} \nonumber \\ &\times \sqrt{\frac{(2\ell''+1)(2\ell'+1)}{8}} \sum_{u,v,w,\tilde{\ell},\tilde{m},\tilde{s},\hat{\ell},\hat{m},\hat{s}} (2 \tilde{\ell}+1) (2 \hat{\ell}+1) (-1)^{\tilde{s}+\hat{s}+\tilde{m}} k_u(\tilde{D}_2) k_v(\tilde{D}_3) k_w(\tilde{D}_1) \nonumber \\ &\qquad\times S(\tilde{\ell},\tilde{D}) \left(\begin{array}{ccc} \ell' & \ell'' & \hat{\ell} \cr -m' & m'' & \hat{m} \end{array}\right) \left(\begin{array}{ccc} \ell' & \ell'' & \hat{\ell} \cr s_{B_1}+s_{B_2}+s_{B_3} & 0 & -\hat{s} \end{array}\right) \nonumber \\ &\qquad\times \left(\begin{array}{ccc} 1 & 1 & \tilde{\ell} \cr -m_u(\tilde{D}_2)&-m_v(\tilde{D}_3)&-\tilde{m} \end{array}\right) \left(\begin{array}{ccc} 1 & 1 & \tilde{\ell} \cr -s_{B_2}&-s_{B_3}&\tilde{s} \end{array}\right) \nonumber \\ &\qquad\times \left(\begin{array}{ccc} 1 & \tilde{\ell} & \hat{\ell} \cr -m_w(\tilde{D}_1)&\tilde{m}&-\hat{m} \end{array}\right) \left(\begin{array}{ccc} 1 & \tilde{\ell} & \hat{\ell} \cr -s_{B_1}&-\tilde{s}&\hat{s} \end{array}\right). \end{align}
As in the rank-2 case, \(S(\tilde{\ell},\tilde{D})\) is a symmetry factor that is unity for a tensor with no symmetries and is
\begin{align} S(\tilde{\ell},\tilde{D}) &= (1+(-1)^{\tilde{\ell}})\frac{2-\delta(\tilde{D_2},\tilde{D_3})}{2} \end{align}
for a tensor symmetric on its last two indices. We don't consider other symmetries because we don't find them in the cases we care about.
| TensorStructure | A Tensor_detail::Structure |
| SparseMatrixType | A sparse matrix fillable by SparseMatrixFiller |
| matrix | The sparse matrix to fill |
| ell_max | The maximum ylm ell value |
| void ylm::TensorYlm::fill_filter | ( | gsl::not_null< SparseMatrixType * > | matrix, |
| size_t | ell_max, | ||
| size_t | number_of_ell_modes_to_kill, | ||
| std::optional< size_t > | half_power | ||
| ) |
Fills a sparse matrix that does a TensorYlm filter operation.
Assumes that \(T^{\tilde A}_{\ell' m'}\) is stored in a Tensor<DataVector>. Multiplying (one plus) the resulting sparse matrix by the Tensor<DataVector> is equivalent to evaluating the right-hand side of Eq. \((\ref{eq:Filter})\).
We assume that the independent components of the Tensor<DataVector> are stored contiguously in memory, in order of the storage_index of the Tensor. However, we do allow for a stride. This means that we can point to memory starting with the first element of \(T^{\tilde A}_{\ell' m'}\) and multiply that by the sparse matrix we compute here, and get the filtered result.
The memory layout here is different than in SpEC. In SpEC, each tensor component is stored in separately-allocated memory, so the SpEC equivalent of the fill_filter function fills \(N^2\) sparse matrices, where \(N\) is the number of independent components of the Tensor. The advantage of the SpEC method is that each sparse matrix is smaller, so sorting elements into the correct order while constructing each sparse matrix is faster (sorting is > linear in the number of matrix elements). The disadvantage of the SpEC method is that evaluating the filter for a single Tensor involves \(N^2\) matrix-vector multiplications, whereas here evaluating the filter involves only one matrix-vector multiplication, which should have more efficient memory access. It is not clear which method is faster overall without more profiling.
The following formulas come from Klinger and Scheel, in prep.
For rank-0 tensors, the expression is simple because there is no change of basis, only a filter based on \(\ell\).
\begin{align} F_{l m \tilde{D}}^{\ell'' m''\tilde{A}} &= \delta(\tilde{D},\tilde{A})\delta_{\ell \ell''}\delta_{m m''} \left[1- \delta(\ell_{\mathrm{cut}}^-\leq \ell \leq \ell_{\mathrm{max}}) g(\ell) \right]. \end{align}
For rank-1 tensors, the expression for \(F_{l m \tilde{D}}^{\ell'' m''\tilde{A}}\) is
\begin{align} F_{l m \tilde{D}}^{\ell'' m''\tilde{A}} &= \delta(\tilde{D},\tilde{A})\delta_{\ell \ell''}\delta_{m m''} \nonumber \\ &- (-1)^{m+m''} (-1)^{\delta(\tilde{D},\mathbf{e}_y)} \delta(\ell,\ell'') \sum_{\ell'=\ell_{\mathrm{cut}}^-}^{\ell_{\mathrm{max}}+1} \frac{2\ell'+1}{2} g(\ell') \nonumber \\ &\times \sum_{j,p,m'} k_j(\tilde{D})k_p(\tilde{A}) \left(\begin{array}{rrr} \ell&\ell'&1\cr -m''&m'&m_j(\tilde{D}) \end{array}\right) \left(\begin{array}{rrr} \ell&\ell'&1\cr -m&m'&m_p(\tilde{A}) \end{array}\right), \end{align}
where the 6-element "matrices" in parentheses are Wigner 3-J symbols, and where
\begin{align} g(\ell') &= \left\{\begin{array}{lr} 1-f(\ell') & \ell' \leq \ell_{\mathrm{cut}}^+,\\ 1 & \ell' > \ell_{\mathrm{cut}}^+. \end{array}\right. \end{align}
For second-rank tensors, the expression for \(F_{l m \tilde{D}}^{\ell'' m''\tilde{A}}\) is
\begin{align} F_{l m \tilde{D}}^{\ell'' m''\tilde{A}} &= \delta(\tilde{D},\tilde{A})\delta_{\ell \ell''}\delta_{m m''} \nonumber \\ &- \frac{1}{4} (-1)^{\delta(\tilde{D}_1,\mathbf{e}_y)} (-1)^{\delta(\tilde{D}_2,\mathbf{e}_y)} \delta_{\ell \ell''} \nonumber \\ &\times \sum_{\ell'=\ell_{\mathrm{cut}}^-}^{\ell_{\mathrm{max}}+2} (2\ell'+1) g(\ell') \sum_{u,v,p,q,\bar{\ell},\tilde{m},\bar{m},m'} (2 \bar{\ell}+1) S(\bar{\ell},\tilde{D}) k_u(\tilde{D}_1) k_v(\tilde{D}_2) k_p(\tilde{A}_1) k_q(\tilde{A}_2) \nonumber \\ &\times \left(\begin{array}{rrr} \ell&\ell'&\bar{\ell}\cr -m&m'&\bar{m} \end{array}\right) \left(\begin{array}{rrr} 1&1&\bar{\ell}\cr m_p(\tilde{A}_1)&m_q(\tilde{A}_2)&-\bar{m} \end{array}\right) \nonumber \\ &\times \left(\begin{array}{rrr} \ell'&\ell&\bar{\ell}\cr -m'&m''&\tilde{m} \end{array}\right) \left(\begin{array}{rrr} 1&1&\bar{\ell}\cr m_u(\tilde{D}_1)&m_v(\tilde{D}_2)&\tilde{m} \end{array}\right), \label{eq:RankTwoTransformWithCut} \end{align}
where the factor \(S(\bar{\ell},\tilde{D})\) is a symmetry factor. For a 2nd-rank tensor with no symmetries, \(S(\bar{\ell},\tilde{D})\) is unity, but for a tensor symmetric in \((\tilde{D_1},\tilde{D_2})\) the matrix elements \(F_{l m \tilde{D}}^{\ell'' m''\tilde{A}}\) are multiplied only by tensor components with \(\tilde{D_1}\geq\tilde{D_2}\) and the symmetry is accounted for by setting
\begin{align} S(\bar{\ell},\tilde{D}) &= (1+(-1)^{\bar{\ell}})\frac{2-\delta(\tilde{D_1},\tilde{D_2})}{2}. \end{align}
For third-rank tensors, the expression for \(F_{l m \tilde{D}}^{\ell'' m''\tilde{A}}\) is
\begin{align} F_{l m \tilde{D}}^{\ell'' m''\tilde{A}} &= \delta(\tilde{D},\tilde{A})\delta_{\ell \ell''}\delta_{m m''} \nonumber \\ &- \frac{1}{8} (-1)^{\delta(\tilde{D}_1,\mathbf{e}_y)} (-1)^{\delta(\tilde{D}_2,\mathbf{e}_y)} (-1)^{\delta(\tilde{D}_3,\mathbf{e}_y)} \delta_{\ell \ell''} \nonumber \\ &\times \sum_{\ell'=\ell_{\mathrm{cut}}^-}^{\ell_{\mathrm{max}}+3} (2\ell'+1) g(\ell') \sum_{u,v,w,p,q,r,\tilde{m},\bar{\ell},\bar{m}, \hat{\ell},\hat{m},\check{m},m'} (2 \bar{\ell}+1) S(\bar{\ell},\tilde{D}) (2 \hat{\ell}+1) \nonumber \\ &\qquad\times k_u(\tilde{D}_2) k_v(\tilde{D}_3) k_w(\tilde{D}_1) k_p(\tilde{A}_2) k_q(\tilde{A}_3) k_r(\tilde{A}_1) (-1)^{\tilde{m}-\bar{m}} \nonumber \\ &\qquad\times \left(\begin{array}{rrr} 1&1&\bar{\ell}\cr m_p(\tilde{A}_2)&m_q(\tilde{A}_3)&-\bar{m} \end{array}\right) \left(\begin{array}{rrr} 1&1&\bar{\ell}\cr m_u(\tilde{D}_2)&m_v(\tilde{D}_3)&\tilde{m} \end{array}\right) \nonumber \\ &\qquad\times \left(\begin{array}{rrr} \ell&\ell'&\hat{\ell}\cr -m&m'&\check{m} \end{array}\right) \left(\begin{array}{rrr} \ell'&\ell&\hat{\ell}\cr -m'&m''&\hat{m} \end{array}\right) \nonumber \\ &\qquad\times \left(\begin{array}{rrr} 1&\bar{\ell}&\hat{\ell}\cr m_r(\tilde{A}_1)&\bar{m}&-\check{m} \end{array}\right) \left(\begin{array}{rrr} 1&\bar{\ell}&\hat{\ell}\cr m_w(\tilde{D}_1)&-\tilde{m}&\hat{m} \end{array}\right). \label{eq:RankThreeTransformWithCut} \end{align}
As in the rank-2 case, \(S(\bar{\ell},\tilde{D})\) is a symmetry factor that is unity for a tensor with no symmetries and is
\begin{align} S(\bar{\ell},\tilde{D}) &= (1+(-1)^{\bar{\ell}})\frac{2-\delta(\tilde{D_2},\tilde{D_3})}{2} \end{align}
for a tensor symmetric on its last two indices. We don't consider other symmetries because we don't find them in the cases we care about.
| TensorStructure | A Tensor_detail::Structure |
| SparseMatrixType | A sparse matrix fillable by SparseMatrixFiller |
| matrix | The sparse matrix to fill |
| ell_max | The maximum ylm ell value. |
| number_of_ell_modes_to_kill | How many top ell modes to set to zero. |
| half_power | The half power \(\sigma\) for more complicated filtering. |
If half_power is std::nullopt, implements a Heaviside filter. Otherwise, the filter is the more complicated one described in the TensorYlm namespace documentation, with \(\sigma\) equal to half_power and \(\ell_{\mathrm{cut}}^+\) equal to \(\ell_{\rm max}\) minus number_of_ell_modes_to_kill.
| void ylm::TensorYlm::fill_sphere_to_cart | ( | gsl::not_null< SparseMatrixType * > | matrix, |
| size_t | ell_max | ||
| ) |
Fills a sparse matrix that does a TensorYlm Cartesian to Spherical operation.
Assumes that the input, \(T^B_{\ell m}\), is stored in a Tensor<DataVector>. Multiplying the resulting sparse matrix by the Tensor<DataVector> is equivalent to evaluating the right-hand side of Eq. \((\ref{eq:S2C})\).
We assume that the independent components of the Tensor<DataVector> are stored contiguously in memory, in order of the storage_index of the Tensor. However, we do allow for a stride. This means that we can point to memory starting with the first element of \(T^B_{\ell m}\) and multiply that by the sparse matrix we compute here, and get the result.
The memory layout here is different than in SpEC. In SpEC, each tensor component is stored in separately-allocated memory, so the SpEC equivalent of the fill_sphere_to_cart function fills \(N^2\) sparse matrices, where \(N\) is the number of independent components of the Tensor. The advantage of the SpEC method is that each sparse matrix is smaller, so sorting elements into the correct order while constructing each sparse matrix is faster (sorting is > linear in the number of matrix elements). The disadvantage of the SpEC method is that evaluating the coefficients for a single Tensor involves \(N^2\) matrix-vector multiplications, whereas here evaluating the coefficients involves only one matrix-vector multiplication, which should have more efficient memory access. It is not clear which method is faster overall without more profiling.
The following formulas come from Klinger and Scheel, in prep.
For rank-1 tensors, the expression for \(C^{\ell' m' \tilde{A}}_{\ell m B}\) is
\begin{align} C^{\ell' m' \tilde{A}}_{\ell m B} &= (-1)^{\delta(s_B,-1)-m}\sqrt{\frac{(2 \ell+1)(2 \ell' + 1)}{2}} \sum_j k_j(\tilde{A}) \left(\begin{array}{ccc} \ell & \ell' & 1 \cr 0 & -s_B & s_B \end{array}\right) \left(\begin{array}{ccc} \ell & \ell' & 1 \cr -m & m' & m_j(\tilde{A}) \end{array}\right), \end{align}
where the 6-element "matrices" in parentheses are Wigner 3-J symbols.
For second-rank tensors, the expression for \(C^{\ell' m' \tilde{A}}_{\ell m B}\) is
\begin{align} C^{\ell' m' \tilde{A}}_{\ell m B} &= \frac{1}{2}(-1)^{\delta(s_{B_1},-1)}(-1)^{\delta(s_{B_2},-1)}(-1)^{m'} \sqrt{(2\ell+1)(2\ell'+1)} \nonumber \\ &\times \sum_{p,q,\bar{\ell},\bar{m},\bar{s}} (2 \bar{\ell}+1) (-1)^{\bar{s}} k_p(\tilde{A}_1) k_q(\tilde{A}_2) S(\bar{\ell},B) \left(\begin{array}{ccc} \ell & \ell' & \bar{\ell} \cr 0&-s_{B_1}-s_{B_2}&-\bar{s} \end{array}\right) \nonumber \\ &\qquad\times \left(\begin{array}{ccc} \ell & \ell' & \bar{\ell} \cr -m & m' & \bar{m} \end{array}\right) \left(\begin{array}{ccc} 1 & 1 & \bar{\ell} \cr m_p(\tilde{A}_1)&m_q(\tilde{A}_2)&-\bar{m} \end{array}\right) \left(\begin{array}{ccc} 1 & 1 & \bar{\ell} \cr s_{B_1}&s_{B_2}&\bar{s} \end{array}\right), \end{align}
where the factor \(S(\bar{\ell},B)\) is a symmetry factor. For a 2nd-rank tensor with no symmetries, \(S(\bar{\ell},B)\) is unity, but for a tensor symmetric in \((B_1,B_2)\) the matrix elements \(C^{\ell' m' \tilde{A}}_{\ell m B}\) are multiplied only by tensor components with \(B_1\geq B_2\) and the symmetry is accounted for by setting
\begin{align} S(\bar{\ell},B) &= (1+(-1)^{\bar{\ell}})\frac{2-\delta(B_1,B_2)}{2}. \end{align}
For third-rank tensors, the expression for \(C^{\ell' m' \tilde{A}}_{\ell m B}\) is
\begin{align} C^{\ell' m' \tilde{A}}_{\ell m B} &= (-1)^{m'} (-1)^{\delta(s_{B_1},-1)} (-1)^{\delta(s_{B_2},-1)} (-1)^{\delta(s_{B_3},-1)} \sqrt{\frac{(2\ell+1)(2\ell'+1)}{8}} \nonumber \\ &\times \sum_{p,q,r,\bar{\ell},\bar{m},\bar{s},\check{\ell},\check{m},\check{s}} (2 \bar{\ell}+1) (2 \check{\ell}+1) (-1)^{\bar{s}+\check{s}-\bar{m}} k_p(\tilde{A}_2) k_q(\tilde{A}_3) k_r(\tilde{A}_1) S(\bar{\ell},B) \nonumber \\ &\qquad\times \left(\begin{array}{ccc} \ell & \ell' & \check{\ell} \cr -m & m' & \check{m} \end{array}\right) \left(\begin{array}{ccc} \ell & \ell' & \check{\ell} \cr 0&-s_{B_1}-s_{B_2}-s_{B_3}&-\check{s} \end{array}\right) \nonumber \\ &\qquad\times \left(\begin{array}{ccc} 1 & 1 & \bar{\ell} \cr m_p(\tilde{A}_2)&m_q(\tilde{A}_3)&-\bar{m} \end{array}\right) \left(\begin{array}{ccc} 1 & 1 & \bar{\ell} \cr s_{B_2}&s_{B_3}&\bar{s} \end{array}\right) \nonumber \\ &\qquad\times \left(\begin{array}{ccc} 1 & \bar{\ell} & \check{\ell} \cr m_r(\tilde{A}_1)&\bar{m}&-\check{m} \end{array}\right) \left(\begin{array}{ccc} 1 & \bar{\ell} & \check{\ell} \cr s_{B_1}&-\bar{s}&\check{s} \end{array}\right), \end{align}
As in the rank-2 case, \(S(\bar{\ell},B)\) is a symmetry factor that is unity for a tensor with no symmetries and is
\begin{align} S(\bar{\ell},B) &= (1+(-1)^{\bar{\ell}})\frac{2-\delta(B_2,B_3)}{2}. \end{align}
for a tensor symmetric on its last two indices. We don't consider other symmetries because we don't find them in the cases we care about.
| TensorStructure | A Tensor_detail::Structure |
| SparseMatrixType | A sparse matrix fillable by SparseMatrixFiller |
| matrix | The sparse matrix to fill |
| ell_max | The maximum ylm ell value |