Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines expressions that represent tensors
6 :
7 : #pragma once
8 :
9 : #include <array>
10 : #include <cassert>
11 : #include <cstddef>
12 : #include <type_traits>
13 : #include <utility>
14 :
15 : #include "DataStructures/Tensor/Expressions/DataTypeSupport.hpp"
16 : #include "DataStructures/Tensor/Expressions/SpatialSpacetimeIndex.hpp"
17 : #include "DataStructures/Tensor/Expressions/TensorExpression.hpp"
18 : #include "Utilities/Algorithm.hpp"
19 : #include "Utilities/ContainerHelpers.hpp"
20 : #include "Utilities/ErrorHandling/Assert.hpp"
21 : #include "Utilities/ForceInline.hpp"
22 : #include "Utilities/Gsl.hpp"
23 : #include "Utilities/Requires.hpp"
24 : #include "Utilities/TMPL.hpp"
25 :
26 : /// \cond
27 : template <typename X, typename Symm, typename IndexList>
28 : class Tensor;
29 : /// \endcond
30 :
31 : namespace tenex {
32 : namespace detail {
33 : /// @{
34 : /// \brief Given a tensor symmetry and the positions of indices where a generic
35 : /// spatial index is used for a spacetime index and the positions where a
36 : /// concrete time index is used, this returns the symmetry after making those
37 : /// indices nonsymmetric with others that do not do the same
38 : ///
39 : /// \details
40 : /// Example: Let `symmetry` be `{1, 2, 2, 2, 1}`. Let
41 : /// `spatial_spacetime_index_positions` be `{1, 2}`: the 2nd and 3rd indices
42 : /// of the tensor are spacetime indices where a generic spatial index is used.
43 : /// Let `time_index_positions` be `{4}`: the 5th (last) index is a spacetime
44 : /// index where a concrete time index is used. The resulting symmetry will
45 : /// reflect that the 2nd and 3rd indices are still symmetric with each other,
46 : /// but no longer symmetric with the 4th index, and the 5th index will no
47 : /// longer be symmetric with the 1st index. Therefore, the returned symmetry
48 : /// will be equivalent to the form of `{4, 3, 3, 2, 1}`.
49 : ///
50 : /// Note: the symmetry returned by this function is not necessarily in the
51 : /// canonical form specified by ::Symmetry. For example, the actual array
52 : /// returned in the example above is `{1, 4, 4, 2, 5}`. As such, this array
53 : /// still encodes which indices are symmetric, but it is not in the canonical
54 : /// form specified by ::Symmetry, which is `{4, 3, 3, 2, 1}`.
55 : ///
56 : /// \param symmetry the input tensor symmetry to transform
57 : /// \param spatial_spacetime_index_positions the positions of the indices of the
58 : /// tensor where a generic spatial index is used for a spacetime index
59 : /// \param time_index_positions the positions of the indices of the tensor where
60 : /// a concrete time index is used for a spacetime index
61 : /// \return the symmetry after making the `spatial_spacetime_index_positions`
62 : /// and `time_index_positions` of `symmetry` nonsymmetric with indices at other
63 : /// positions
64 : template <
65 : size_t NumIndices, size_t NumSpatialSpacetimeIndices,
66 : size_t NumConcreteTimeIndices,
67 : Requires<(NumIndices >= 2 and (NumSpatialSpacetimeIndices != 0 or
68 : NumConcreteTimeIndices != 0))> = nullptr>
69 : SPECTRE_ALWAYS_INLINE constexpr std::array<std::int32_t, NumIndices>
70 : get_transformed_spacetime_symmetry(
71 : const std::array<std::int32_t, NumIndices>& symmetry,
72 : const std::array<size_t, NumSpatialSpacetimeIndices>&
73 : spatial_spacetime_index_positions,
74 : const std::array<size_t, NumConcreteTimeIndices>& time_index_positions) {
75 : std::array<std::int32_t, NumIndices> transformed_spacetime_symmetry{};
76 : const std::int32_t max_symm_value =
77 : static_cast<std::int32_t>(*alg::max_element(symmetry));
78 : for (size_t i = 0; i < NumIndices; i++) {
79 : gsl::at(transformed_spacetime_symmetry, i) = gsl::at(symmetry, i);
80 : }
81 : for (size_t i = 0; i < NumSpatialSpacetimeIndices; i++) {
82 : gsl::at(transformed_spacetime_symmetry,
83 : gsl::at(spatial_spacetime_index_positions, i)) += max_symm_value;
84 : }
85 : for (size_t i = 0; i < NumConcreteTimeIndices; i++) {
86 : gsl::at(transformed_spacetime_symmetry, gsl::at(time_index_positions, i)) +=
87 : 2 * max_symm_value;
88 : }
89 :
90 : return transformed_spacetime_symmetry;
91 : }
92 :
93 : template <size_t NumIndices, size_t NumSpatialSpacetimeIndices,
94 : size_t NumConcreteTimeIndices,
95 : Requires<(NumIndices < 2 or (NumSpatialSpacetimeIndices == 0 and
96 : NumConcreteTimeIndices == 0))> = nullptr>
97 : SPECTRE_ALWAYS_INLINE constexpr std::array<std::int32_t, NumIndices>
98 : get_transformed_spacetime_symmetry(
99 : const std::array<std::int32_t, NumIndices>& symmetry,
100 : const std::array<size_t, NumSpatialSpacetimeIndices>&
101 : /*spatial_spacetime_index_positions*/,
102 : const std::array<size_t, NumConcreteTimeIndices>&
103 : /*time_index_positions*/) {
104 : return symmetry;
105 : }
106 : /// @}
107 :
108 : /// \ingroup TensorExpressionsGroup
109 : /// \brief Helper struct for computing the new canonical symmetry of a tensor
110 : /// after generic spatial indices are used for any of the tensor's spacetime
111 : /// indices
112 : ///
113 : /// \details This is relevant in cases where a tensor has spacetime indices that
114 : /// are symmetric but generic spatial indices are used for a non-empty subset of
115 : /// those symmetric spacetime indices. For example, if we have some rank 3
116 : /// tensor with the first index being spatial and the 2nd and third indices
117 : /// spacetime and symmetric, but a generic spatial index is used for the 2nd
118 : /// index, the "result" of the single tensor expression, \f$R_{ija}\f$, is a
119 : /// rank 3 tensor whose 2nd and 3rd indices are no longer symmetric.
120 :
121 : /// Given that some `Tensor` named `R` that represents the tensor in the above
122 : /// example, the symmetry of the `Tensor` is `[2, 1, 1]`, but the computed
123 : /// symmetry of the `TensorAsExpression` that represents it will have symmetry
124 : /// `[3, 2, 1]` to reflect this loss of symmetry.
125 : ///
126 : /// Evaluating the "result" symmetry here in `TensorAsExpression`, at the leaves
127 : /// of the expression tree, enables the propagation of this symmetry up the tree
128 : /// to the other expression types. By determining each tensor's "result"
129 : /// symmetry at the leaves, the expressions at internal nodes of the tree can
130 : /// have their individual symmetries determined without having to each consider
131 : /// whether their operand(s) are expression(s) that have spacetime indices where
132 : /// generic spatial indices were used.
133 : ///
134 : /// \tparam SymmList the ::Symmetry of the Tensor represented by the expression
135 : /// \tparam TensorIndexTypeList the \ref SpacetimeIndex "TensorIndexType"'s of
136 : /// the Tensor represented by the expression
137 : /// \tparam TensorIndexList the generic indices of the Tensor represented by the
138 : /// expression
139 : template <typename SymmList, typename TensorIndexTypeList,
140 : typename TensorIndexList,
141 : size_t NumIndices = tmpl::size<SymmList>::value,
142 : typename IndexSequence = std::make_index_sequence<NumIndices>>
143 : struct TensorAsExpressionSymm;
144 :
145 : template <template <typename...> class SymmList, typename... Symm,
146 : typename TensorIndexTypeList, typename TensorIndexList,
147 : size_t NumIndices, size_t... Ints>
148 : struct TensorAsExpressionSymm<SymmList<Symm...>, TensorIndexTypeList,
149 : TensorIndexList, NumIndices,
150 : std::index_sequence<Ints...>> {
151 : static constexpr auto symm = get_transformed_spacetime_symmetry<NumIndices>(
152 : {{Symm::value...}},
153 : get_spatial_spacetime_index_positions<TensorIndexTypeList,
154 : TensorIndexList>(),
155 : get_time_index_positions<TensorIndexList>());
156 : using type = Symmetry<symm[Ints]...>;
157 : };
158 : } // namespace detail
159 :
160 : /// \ingroup TensorExpressionsGroup
161 : /// \brief Defines an expression representing a Tensor
162 : ///
163 : /// \details
164 : /// In order to represent a tensor as an expression, instead of having Tensor
165 : /// derive off of TensorExpression, a TensorAsExpression derives off of
166 : /// TensorExpression and contains a pointer to a Tensor. The reason having
167 : /// Tensor derive off of TensorExpression is problematic is that the index
168 : /// structure is part of the type of the TensorExpression, so every possible
169 : /// permutation and combination of indices must be derived from. For a rank 3
170 : /// tensor, this is already over 500 base classes, which the Intel compiler
171 : /// takes too long to compile.
172 : ///
173 : /// For details on aliases and members defined in this class, as well as general
174 : /// `TensorExpression` terminology used in its members' documentation, see
175 : /// documentation for `TensorExpression`.
176 : ///
177 : /// \tparam T the type of Tensor being represented as an expression
178 : /// \tparam ArgsList the tensor indices, e.g. `_a` and `_b` in `F(_a, _b)`
179 : template <typename T, typename ArgsList>
180 1 : struct TensorAsExpression;
181 :
182 : template <typename X, template <typename...> class Symm, typename... SymmValues,
183 : template <typename...> class IndexList, typename... Indices,
184 : template <typename...> class ArgsList, typename... Args>
185 0 : struct TensorAsExpression<Tensor<X, Symm<SymmValues...>, IndexList<Indices...>>,
186 : ArgsList<Args...>>
187 : : public TensorExpression<TensorAsExpression<Tensor<X, Symm<SymmValues...>,
188 : IndexList<Indices...>>,
189 : ArgsList<Args...>>,
190 : X,
191 : typename detail::TensorAsExpressionSymm<
192 : Symm<SymmValues...>, IndexList<Indices...>,
193 : ArgsList<Args...>>::type,
194 : IndexList<Indices...>, ArgsList<Args...>> {
195 : static_assert(detail::is_supported_tensor_datatype_v<X>,
196 : "TensorExpressions currently only support Tensors whose data "
197 : "type is double, std::complex<double>, DataVector, or "
198 : "ComplexDataVector. It is possible to add support for other "
199 : "data types that are supported by Tensor.");
200 : // `Symmetry` currently prevents this because antisymmetries are not currently
201 : // supported for `Tensor`s. This check is repeated here because if
202 : // antisymmetries are later supported for `Tensor`, using antisymmetries in
203 : // `TensorExpression`s will not automatically work. The implementations of the
204 : // derived `TensorExpression` types assume no antisymmetries (assume positive
205 : // `Symmetry` values), so support for antisymmetries in `TensorExpression`s
206 : // will still need to be implemented.
207 : static_assert((... and (SymmValues::value > 0)),
208 : "Anti-symmetric Tensors are currently not supported by "
209 : "TensorExpressions.");
210 :
211 : // === Index properties ===
212 : /// The type of the data being stored in the result of the expression
213 1 : using type = X;
214 : /// The list of \ref SpacetimeIndex "TensorIndexType"s of the result of the
215 : /// expression
216 1 : using symmetry = typename detail::TensorAsExpressionSymm<
217 : Symm<SymmValues...>, IndexList<Indices...>, ArgsList<Args...>>::type;
218 : /// The list of \ref SpacetimeIndex "TensorIndexType"s of the result of the
219 : /// expression
220 1 : using index_list = IndexList<Indices...>;
221 : /// The list of generic `TensorIndex`s of the result of the expression
222 1 : using args_list = ArgsList<Args...>;
223 : /// The number of tensor indices in the result of the expression
224 1 : static constexpr auto num_tensor_indices = tmpl::size<index_list>::value;
225 :
226 : // === Expression subtree properties ===
227 : /// The number of arithmetic tensor operations done in the subtree for the
228 : /// left operand, which is 0 because this is a leaf expression
229 1 : static constexpr size_t num_ops_left_child = 0;
230 : /// The number of arithmetic tensor operations done in the subtree for the
231 : /// right operand, which is 0 because this is a leaf expression
232 1 : static constexpr size_t num_ops_right_child = 0;
233 : /// The total number of arithmetic tensor operations done in this expression's
234 : /// whole subtree, which is 0 because this is a leaf expression
235 1 : static constexpr size_t num_ops_subtree = 0;
236 : /// The height of this expression's node in the expression tree relative to
237 : /// the closest `TensorAsExpression` leaf in its subtree. Because this
238 : /// expression type is that leaf, the height for this type will always be 0.
239 1 : static constexpr size_t height_relative_to_closest_tensor_leaf_in_subtree = 0;
240 :
241 : // === Properties for splitting up subexpressions along the primary path ===
242 : // These definitions only have meaning if this expression actually ends up
243 : // being along the primary path that is taken when evaluating the whole tree.
244 : // See documentation for `TensorExpression` for more details.
245 : /// If on the primary path, whether or not the expression is an ending point
246 : /// of a leg
247 1 : static constexpr bool is_primary_end = true;
248 : /// If on the primary path, this is the remaining number of arithmetic tensor
249 : /// operations that need to be done in the subtree of the child along the
250 : /// primary path, given that we will have already computed the whole subtree
251 : /// at the next lowest leg's starting point. This is just 0 because this
252 : /// expression is a leaf.
253 1 : static constexpr size_t num_ops_to_evaluate_primary_left_child = 0;
254 : /// If on the primary path, this is the remaining number of arithmetic tensor
255 : /// operations that need to be done in the right operand's subtree. This is
256 : /// just 0 because this expression is a leaf.
257 1 : static constexpr size_t num_ops_to_evaluate_primary_right_child = 0;
258 : /// If on the primary path, this is the remaining number of arithmetic tensor
259 : /// operations that need to be done for this expression's subtree, given that
260 : /// we will have already computed the subtree at the next lowest leg's
261 : /// starting point. This is just 0 because this expression is a leaf.
262 1 : static constexpr size_t num_ops_to_evaluate_primary_subtree = 0;
263 : /// If on the primary path, whether or not the expression is a starting point
264 : /// of a leg
265 1 : static constexpr bool is_primary_start = false;
266 : /// If on the primary path, whether or not the expression's child along the
267 : /// primary path is a subtree that contains a starting point of a leg along
268 : /// the primary path. This is always falls because this expression is a leaf.
269 1 : static constexpr bool primary_child_subtree_contains_primary_start = false;
270 : /// If on the primary path, whether or not this subtree contains a starting
271 : /// point of a leg along the primary path
272 1 : static constexpr bool primary_subtree_contains_primary_start =
273 : is_primary_start;
274 :
275 : /// Construct an expression from a Tensor
276 1 : explicit TensorAsExpression(
277 : const Tensor<X, Symm<SymmValues...>, IndexList<Indices...>>& t)
278 : : t_(&t) {}
279 0 : ~TensorAsExpression() override = default;
280 :
281 : /// \brief Assert that the LHS tensor of the equation is not equal to the
282 : /// `Tensor` represented by this expression
283 : template <typename LhsTensor>
284 1 : SPECTRE_ALWAYS_INLINE void assert_lhs_tensor_not_in_rhs_expression(
285 : const gsl::not_null<LhsTensor*> lhs_tensor) const {
286 : (void)lhs_tensor;
287 : ASSERT(static_cast<const void*>(&(*t_)) != &(*lhs_tensor),
288 : "The LHS Tensor cannot also be in the RHS expression in the call to "
289 : "tenex::evaluate(). Use tenex::update() instead.");
290 : }
291 :
292 : /// \brief If the LHS tensor is the tensor represented by this expression,
293 : /// assert that the order of the generic indices are the same
294 : ///
295 : /// \tparam LhsTensorIndices the list of generic `TensorIndex`s of the LHS
296 : /// result `Tensor` being computed
297 : /// \param lhs_tensor the LHS result `Tensor` being computed
298 : template <typename LhsTensorIndices, typename LhsTensor>
299 1 : SPECTRE_ALWAYS_INLINE void assert_lhs_tensorindices_same_in_rhs(
300 : const gsl::not_null<LhsTensor*> lhs_tensor) const {
301 : if (static_cast<const void*>(&(*t_)) == &(*lhs_tensor)) {
302 : ASSERT(
303 : (std::is_same_v<LhsTensorIndices, args_list>),
304 : "The LHS Tensor was also found in the RHS tensor expression, but the "
305 : "generic index order used for it on the RHS does not match the order "
306 : "used for it on the LHS.");
307 : }
308 : }
309 :
310 : /// \brief Get the size of a component from the `Tensor` contained by this
311 : /// expression
312 : ///
313 : /// \return the size of a component from the `Tensor` contained by this
314 : /// expression
315 1 : SPECTRE_ALWAYS_INLINE size_t get_rhs_tensor_component_size() const {
316 : return get_size((*t_)[0]);
317 : }
318 :
319 : /// \brief Returns the value of the contained tensor's multi-index
320 : ///
321 : /// \param multi_index the multi-index of the tensor component to retrieve
322 : /// \return the value of the component at `multi_index` in the tensor
323 1 : SPECTRE_ALWAYS_INLINE decltype(auto) get(
324 : const std::array<size_t, num_tensor_indices>& multi_index) const {
325 : return t_->get(multi_index);
326 : }
327 :
328 : /// \brief Returns the value of the contained tensor's multi-index
329 : ///
330 : /// \param multi_index the multi-index of the tensor component to retrieve
331 : /// \return the value of the component at `multi_index` in the tensor
332 : template <typename ResultType>
333 1 : SPECTRE_ALWAYS_INLINE decltype(auto) get_primary(
334 : const ResultType& /*result_component*/,
335 : const std::array<size_t, num_tensor_indices>& multi_index) const {
336 : return t_->get(multi_index);
337 : }
338 :
339 : // This expression is a leaf and therefore will never be the start of a leg
340 : // to evaluate in a split tree, which is enforced by
341 : // `is_primary_start == false`. Therefore, this function should never be
342 : // called on this expression type.
343 : template <typename ResultType>
344 0 : SPECTRE_ALWAYS_INLINE void evaluate_primary_subtree(
345 : ResultType& result_component,
346 : const std::array<size_t, num_tensor_indices>&) const = delete;
347 :
348 : /// Retrieve the i'th entry of the Tensor being held
349 1 : SPECTRE_ALWAYS_INLINE type operator[](const size_t i) const {
350 : return t_->operator[](i);
351 : }
352 :
353 : private:
354 : /// `Tensor` represented by this expression
355 1 : const Tensor<X, Symm<SymmValues...>, IndexList<Indices...>>* t_ = nullptr;
356 : };
357 : } // namespace tenex
|
