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