SpECTRE Documentation Coverage Report
Current view: top level - DataStructures/Tensor/Expressions - AddSubtract.hpp Hit Total Coverage
Commit: c7c91c0aa25aeda7dc9e4c735164d2bb71ebe72e Lines: 46 51 90.2 %
Date: 2024-03-28 06:01:39
Legend: Lines: hit not hit

          Line data    Source code
       1           1 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : /// \file
       5             : /// Defines ET for adding and subtracting tensors
       6             : 
       7             : #pragma once
       8             : 
       9             : #include <array>
      10             : #include <complex>
      11             : #include <cstddef>
      12             : #include <iterator>
      13             : #include <limits>
      14             : #include <type_traits>
      15             : #include <utility>
      16             : 
      17             : #include "DataStructures/DataVector.hpp"
      18             : #include "DataStructures/Tensor/Expressions/DataTypeSupport.hpp"
      19             : #include "DataStructures/Tensor/Expressions/IndexPropertyCheck.hpp"
      20             : #include "DataStructures/Tensor/Expressions/NumberAsExpression.hpp"
      21             : #include "DataStructures/Tensor/Expressions/SpatialSpacetimeIndex.hpp"
      22             : #include "DataStructures/Tensor/Expressions/TensorExpression.hpp"
      23             : #include "DataStructures/Tensor/Expressions/TensorIndex.hpp"
      24             : #include "DataStructures/Tensor/Expressions/TensorIndexTransformation.hpp"
      25             : #include "DataStructures/Tensor/Expressions/TimeIndex.hpp"
      26             : #include "Utilities/Algorithm.hpp"
      27             : #include "Utilities/ForceInline.hpp"
      28             : #include "Utilities/MakeArray.hpp"
      29             : #include "Utilities/Requires.hpp"
      30             : #include "Utilities/TMPL.hpp"
      31             : 
      32             : /// \cond
      33             : namespace tenex {
      34             : template <typename T1, typename T2, typename ArgsList1, typename ArgsList2,
      35             :           int Sign>
      36             : struct AddSub;
      37             : }  // namespace tenex
      38             : template <typename Derived, typename DataType, typename Symm,
      39             :           typename IndexList, typename Args, typename ReducedArgs>
      40             : struct TensorExpression;
      41             : /// \endcond
      42             : 
      43           1 : namespace tenex {
      44             : namespace detail {
      45             : /// \brief Computes the rearranged symmetry of one operand according to the
      46             : /// generic index order of the other operand
      47             : ///
      48             : /// \details
      49             : /// Here is an example of what the algorithm does:
      50             : ///
      51             : /// Given \f$R_{abc} + S_{cab}\f$, reorder \f$S\f$' symmetry according to
      52             : /// \f$R\f$'s index order
      53             : /// `tensorindex_transformation`:
      54             : /// \code
      55             : /// {1, 2, 0} // positions of R's indices {a, b, c} in S' indices {c, a, b}
      56             : /// \endcode
      57             : /// `input_symm`:
      58             : /// \code
      59             : /// {2, 2, 1} // S' symmetry, where c and a are symmetric indices
      60             : /// \endcode
      61             : /// returned equivalent `output_symm`:
      62             : /// \code
      63             : /// {2, 1, 2} // S' symmetry (input_symm) rearranged to R's index order (abc)
      64             : /// \endcode
      65             : ///
      66             : /// One special case scenario to note is when concrete time indices are
      67             : /// involved in the transformation. Consider transforming the symmetry for
      68             : /// some tensor \f$S_{ab}\f$ to the index order of another tensor
      69             : /// \f$R_{btat}\f$. This would be necessary in an expression such as
      70             : /// \f$R_{btat} + S_{ab}\f$. The transformation of the symmetry for \f$S\f$
      71             : /// according to the index order of \f$R\f$ cannot simply be the list of
      72             : /// positions of \f$R\f$'s indices in \f$S\f$' indices, as \f$S\f$ does not
      73             : /// contain all of \f$R\f$'s indices, because it has no time indices. To handle
      74             : /// cases like this, a placeholder value for the position of any time index
      75             : /// must be substituted for an actual position, since one may not exist. In this
      76             : /// example, the proper input transformation (`tensorindex_transformation`)
      77             : /// would need to be `{1, PLACEHOLDER_VALUE, 0, PLACEHOLDER_VALUE}`, where
      78             : /// `PLACEHOLDER_VALUE` is defined by
      79             : /// `TensorIndexTransformation_detail::time_index_position_placeholder`. `1` and
      80             : /// `0` are the positions of \f$b\f$ and \f$a\f$ in \f$S\f$, and the placeholder
      81             : /// is used for the positions of time indices. In computing the output
      82             : /// transformed symmetry, the function will insert a `0` at each position where
      83             : /// this placeholder is found in the transformation. For example, if
      84             : /// `input_symm` is `{2, 1}`, the returned output multi-index will be
      85             : /// `{1, 0, 2, 0}`. Note that the symmetry returned by this function is not
      86             : /// necessarily in the canonical form defined by ::Symmetry. This example with
      87             : /// time indices is an example of this, as `0` is not a permitted ::Symmetry
      88             : /// value, and the canonical form would have increasing symmetry values from
      89             : /// right to left. In addition, even though the time indices in the rearranged
      90             : /// symmetry will have the same symmetry value (`0`), this bears no impact on
      91             : /// `get_addsub_symm`'s computation of the symmetry of the tensor resulting
      92             : /// from the addition or subtraction.
      93             : ///
      94             : /// \tparam NumIndicesIn the number of indices in the operand whose symmetry is
      95             : /// being transformed
      96             : /// \tparam NumIndicesOut the number of indices in the other operand whose index
      97             : /// order is the order the input operand symmetry is being transformed to
      98             : /// \param input_symm the input operand symmetry to transform
      99             : /// \param tensorindex_transformation the positions of the other operand's
     100             : /// generic indices in the generic indices of the operand whose symmetry is
     101             : /// being transformed (see details)
     102             : /// \return the input operand symmetry rearranged according to the generic index
     103             : /// order of the other operand
     104             : template <size_t NumIndicesIn, size_t NumIndicesOut>
     105             : SPECTRE_ALWAYS_INLINE constexpr std::array<std::int32_t, NumIndicesOut>
     106             : transform_addsub_symm(
     107             :     const std::array<std::int32_t, NumIndicesIn>& input_symm,
     108             :     const std::array<size_t, NumIndicesOut>& tensorindex_transformation) {
     109             :   std::array<std::int32_t, NumIndicesOut> output_symm =
     110             :       make_array<NumIndicesOut, std::int32_t>(0);
     111             :   for (size_t i = 0; i < NumIndicesOut; i++) {
     112             :     gsl::at(output_symm, i) =
     113             :         (gsl::at(tensorindex_transformation, i) ==
     114             :          TensorIndexTransformation_detail::time_index_position_placeholder)
     115             :             ? 0
     116             :             : gsl::at(input_symm, gsl::at(tensorindex_transformation, i));
     117             :   }
     118             :   return output_symm;
     119             : }
     120             : 
     121             : /// @{
     122             : /// \ingroup TensorExpressionsGroup
     123             : /// \brief Returns the canonical symmetry of the tensor resulting from
     124             : /// adding or subtracting two tensors, according to their symmetries
     125             : ///
     126             : /// \details The canonical symmetry returned follows the convention defined by
     127             : /// ::Symmetry: symmetry values are in ascending order from right to left. If
     128             : /// the convention implemented by ::Symmetry changes, this function will also
     129             : /// need to be updated to match the new convention. The ::Symmetry metafunction
     130             : /// could instead be used on the result of this function, but that would
     131             : /// introduce avoidable and unnecessary extra computations, so it is not used.
     132             : ///
     133             : /// This function treats the two input symmetries as aligned (i.e. each position
     134             : /// of `symm1` and `symm2` corresponds to a shared generic index at that
     135             : /// position). The resultant symmetry is determined as follows: indices that are
     136             : /// symmetric in both input symmetries are also symmetric in the resultant
     137             : /// tensor.
     138             : ///
     139             : /// \param symm1 the symmetry of the first tensor being added or subtracted
     140             : /// \param symm2 the symmetry of the second tensor being added or subtracted
     141             : /// \return the canonical symmetry of the tensor resulting from adding or
     142             : /// subtracting two tensors
     143             : template <size_t NumIndices, Requires<(NumIndices >= 2)> = nullptr>
     144             : constexpr std::array<std::int32_t, NumIndices> get_addsub_symm(
     145             :     const std::array<std::int32_t, NumIndices>& symm1,
     146             :     const std::array<std::int32_t, NumIndices>& symm2) {
     147             :   constexpr std::int32_t max_int = std::numeric_limits<std::int32_t>::max();
     148             :   std::array<std::int32_t, NumIndices> addsub_symm =
     149             :       make_array<NumIndices>(max_int);
     150             :   size_t right_index = NumIndices - 1;
     151             :   std::int32_t symm_value_to_set = 1;
     152             : 
     153             :   while (right_index < NumIndices) {
     154             :     std::int32_t symm1_value_to_find = symm1[right_index];
     155             :     std::int32_t symm2_value_to_find = symm2[right_index];
     156             :     // if we haven't yet set right_index for the resultant symmetry
     157             :     if (addsub_symm[right_index] == max_int) {
     158             :       addsub_symm[right_index] = symm_value_to_set;
     159             :       for (size_t left_index = right_index - 1; left_index < NumIndices;
     160             :            left_index--) {
     161             :         // if left_index of the resultant symmetry is not yet set and we've
     162             :         // found a common symmetry between symm1 and symm2 at this index
     163             :         if (addsub_symm[left_index] == max_int and
     164             :             symm1[left_index] == symm1_value_to_find and
     165             :             symm2[left_index] == symm2_value_to_find) {
     166             :           addsub_symm[left_index] = symm_value_to_set;
     167             :         }
     168             :       }
     169             :       symm_value_to_set++;
     170             :     }
     171             :     right_index--;
     172             :   }
     173             : 
     174             :   return addsub_symm;
     175             : }
     176             : 
     177             : template <size_t NumIndices, Requires<(NumIndices == 1)> = nullptr>
     178             : constexpr std::array<std::int32_t, NumIndices> get_addsub_symm(
     179             :     const std::array<std::int32_t, NumIndices>& /*symm1*/,
     180             :     const std::array<std::int32_t, NumIndices>& /*symm2*/) {
     181             :   // return {{1}} instead of symm1 in case symm1 is not in the canonical form
     182             :   return {{1}};
     183             : }
     184             : 
     185             : template <size_t NumIndices, Requires<(NumIndices == 0)> = nullptr>
     186             : constexpr std::array<std::int32_t, NumIndices> get_addsub_symm(
     187             :     const std::array<std::int32_t, NumIndices>& symm1,
     188             :     const std::array<std::int32_t, NumIndices>& /*symm2*/) {
     189             :   return symm1;
     190             : }
     191             : /// @}
     192             : 
     193             : /// \ingroup TensorExpressionsGroup
     194             : /// \brief Helper struct for computing the canonical symmetry of the tensor
     195             : /// resulting from adding or subtracting two tensors, according to their
     196             : /// symmetries and generic index orders
     197             : ///
     198             : /// \details The resultant symmetry (`type`) values correspond to the index
     199             : /// order of the first tensor operand being added or subtracted:
     200             : /// `TensorIndexList1`.
     201             : ///
     202             : /// \tparam SymmList1 the ::Symmetry of the first operand
     203             : /// \tparam SymmList2 the ::Symmetry of the second operand
     204             : /// \tparam TensorIndexList1 the generic indices of the first operand
     205             : /// \tparam TensorIndexList2 the generic indices of the second operand
     206             : template <typename SymmList1, typename SymmList2, typename TensorIndexList1,
     207             :           typename TensorIndexList2,
     208             :           size_t NumIndices1 = tmpl::size<SymmList1>::value,
     209             :           size_t NumIndices2 = tmpl::size<SymmList2>::value,
     210             :           typename IndexSequence1 = std::make_index_sequence<NumIndices1>>
     211             : struct AddSubSymmetry;
     212             : 
     213             : template <
     214             :     template <typename...> class SymmList1, typename... Symm1,
     215             :     template <typename...> class SymmList2, typename... Symm2,
     216             :     template <typename...> class TensorIndexList1, typename... TensorIndices1,
     217             :     template <typename...> class TensorIndexList2, typename... TensorIndices2,
     218             :     size_t NumIndices1, size_t NumIndices2, size_t... Ints1>
     219             : struct AddSubSymmetry<SymmList1<Symm1...>, SymmList2<Symm2...>,
     220             :                       TensorIndexList1<TensorIndices1...>,
     221             :                       TensorIndexList2<TensorIndices2...>, NumIndices1,
     222             :                       NumIndices2, std::index_sequence<Ints1...>> {
     223             :   static constexpr std::array<size_t, NumIndices1> tensorindex_values1 = {
     224             :       {TensorIndices1::value...}};
     225             :   static constexpr std::array<size_t, NumIndices2> tensorindex_values2 = {
     226             :       {TensorIndices2::value...}};
     227             :   // positions of tensorindex_values1 in tensorindex_values2
     228             :   static constexpr std::array<size_t, NumIndices1> op2_to_op1_map =
     229             :       ::tenex::compute_tensorindex_transformation(tensorindex_values2,
     230             :                                                   tensorindex_values1);
     231             : 
     232             :   static constexpr std::array<std::int32_t, NumIndices1> symm1 = {
     233             :       {Symm1::value...}};
     234             :   static constexpr std::array<std::int32_t, NumIndices2> symm2 = {
     235             :       {Symm2::value...}};
     236             :   // 2nd argument is symm2 rearranged according to `TensorIndexList1` order
     237             :   // so that the two symmetry arguments to `get_addsub_symm` are aligned
     238             :   // w.r.t. their generic index orders
     239             :   static constexpr std::array<std::int32_t, NumIndices1> addsub_symm =
     240             :       get_addsub_symm(symm1, transform_addsub_symm(symm2, op2_to_op1_map));
     241             : 
     242             :   using type = tmpl::integral_list<std::int32_t, addsub_symm[Ints1]...>;
     243             : };
     244             : 
     245             : /// \ingroup TensorExpressionsGroup
     246             : /// \brief Helper struct for defining the symmetry, index list, and
     247             : /// generic index list of the tensor resulting from adding or
     248             : /// subtracting two tensor expressions
     249             : ///
     250             : /// \tparam T1 the first tensor expression operand
     251             : /// \tparam T2 the second tensor expression operand
     252             : template <typename T1, typename T2>
     253             : struct AddSubType {
     254             :   static_assert(std::is_base_of_v<Expression, T1> and
     255             :                     std::is_base_of_v<Expression, T2>,
     256             :                 "Parameters to AddSubType must be TensorExpressions");
     257             :   using type =
     258             :       typename get_binop_datatype<typename T1::type, typename T2::type>::type;
     259             :   using symmetry =
     260             :       typename AddSubSymmetry<typename T1::symmetry, typename T2::symmetry,
     261             :                                    typename T1::args_list,
     262             :                                    typename T2::args_list>::type;
     263             :   using index_list = typename T1::index_list;
     264             :   using tensorindex_list = typename T1::args_list;
     265             : };
     266             : }  // namespace detail
     267             : 
     268             : /// \ingroup TensorExpressionsGroup
     269             : /// \brief Defines the tensor expression representing the addition or
     270             : /// subtraction of two tensor expressions
     271             : ///
     272             : /// \details
     273             : /// For details on aliases and members defined in this class, as well as general
     274             : /// `TensorExpression` terminology used in its members' documentation, see
     275             : /// documentation for `TensorExpression`.
     276             : ///
     277             : /// \tparam T1 the left operand expression
     278             : /// \tparam T2 the right operand expression
     279             : /// \tparam ArgsList1 generic `TensorIndex`s of the left operand
     280             : /// \tparam ArgsList2 generic `TensorIndex`s of the right operand
     281             : /// \tparam Sign the sign of the operation selected, 1 for addition or -1 for
     282             : /// subtraction
     283             : template <typename T1, typename T2, typename ArgsList1, typename ArgsList2,
     284             :           int Sign>
     285           1 : struct AddSub;
     286             : 
     287             : template <typename T1, typename T2, template <typename...> class ArgsList1,
     288             :           template <typename...> class ArgsList2, typename... Args1,
     289             :           typename... Args2, int Sign>
     290           0 : struct AddSub<T1, T2, ArgsList1<Args1...>, ArgsList2<Args2...>, Sign>
     291             :     : public TensorExpression<
     292             :           AddSub<T1, T2, ArgsList1<Args1...>, ArgsList2<Args2...>, Sign>,
     293             :           typename detail::AddSubType<T1, T2>::type,
     294             :           typename detail::AddSubType<T1, T2>::symmetry,
     295             :           typename detail::AddSubType<T1, T2>::index_list,
     296             :           typename detail::AddSubType<T1, T2>::tensorindex_list> {
     297             :   static_assert(
     298             :       detail::tensorexpression_binop_datatypes_are_supported_v<T1, T2>,
     299             :       "Cannot add or subtract the given TensorExpression types with the given "
     300             :       "data types. This can occur from e.g. trying to add a Tensor with data "
     301             :       "type double and a Tensor with data type DataVector.");
     302             :   static_assert(
     303             :       not((std::is_same_v<T1, NumberAsExpression<std::complex<double>>> and
     304             :            std::is_same_v<typename T2::type, DataVector>) or
     305             :           (std::is_same_v<T2, NumberAsExpression<std::complex<double>>> and
     306             :            std::is_same_v<typename T1::type, DataVector>)),
     307             :       "Cannot perform addition and subtraction between a std::complex<double> "
     308             :       "and a TensorExpression whose data type is DataVector because Blaze does "
     309             :       "not support addition and subtraction between std::complex<double> and "
     310             :       "DataVector.");
     311             :   static_assert(
     312             :       detail::IndexPropertyCheck<typename T1::index_list,
     313             :                                  typename T2::index_list, ArgsList1<Args1...>,
     314             :                                  ArgsList2<Args2...>>::value,
     315             :       "You are attempting to add indices of different types, e.g. T^a_b + "
     316             :       "S^b_a, which doesn't make sense. The indices may also be in different "
     317             :       "frames, different types (spatial vs. spacetime) or of different "
     318             :       "dimension.");
     319             :   static_assert(Sign == 1 or Sign == -1,
     320             :                 "Invalid Sign provided for addition or subtraction of Tensor "
     321             :                 "elements. Sign must be 1 (addition) or -1 (subtraction).");
     322             : 
     323             :   // === Index properties ===
     324             :   /// The type of the data being stored in the result of the expression
     325           1 :   using type = typename detail::AddSubType<T1, T2>::type;
     326             :   /// The ::Symmetry of the result of the expression
     327           1 :   using symmetry = typename detail::AddSubType<T1, T2>::symmetry;
     328             :   /// The list of \ref SpacetimeIndex "TensorIndexType"s of the result of the
     329             :   /// expression
     330           1 :   using index_list = typename detail::AddSubType<T1, T2>::index_list;
     331             :   /// The list of generic `TensorIndex`s of the result of the
     332             :   /// expression
     333           1 :   using args_list = typename T1::args_list;
     334             :   /// The number of tensor indices in the result of the expression. This also
     335             :   /// doubles as the left operand's number of indices.
     336           1 :   static constexpr auto num_tensor_indices = tmpl::size<index_list>::value;
     337             :   /// The number of tensor indices in the right operand expression
     338           1 :   static constexpr auto num_tensor_indices_op2 = sizeof...(Args2);
     339             :   /// Mapping from the left operand's index order to the right operand's index
     340             :   /// order
     341             :   static constexpr std::array<size_t, num_tensor_indices_op2>
     342           1 :       operand_index_transformation =
     343             :           compute_tensorindex_transformation<num_tensor_indices,
     344             :                                              num_tensor_indices_op2>(
     345             :               {{Args1::value...}}, {{Args2::value...}});
     346             :   /// Positions of indices in first operand where generic spatial indices are
     347             :   /// used for spacetime indices
     348           1 :   static constexpr auto op1_spatial_spacetime_index_positions =
     349             :       detail::get_spatial_spacetime_index_positions<typename T1::index_list,
     350             :                                                     ArgsList1<Args1...>>();
     351             :   /// Positions of indices in second operand where generic spatial indices are
     352             :   /// used for spacetime indices
     353           1 :   static constexpr auto op2_spatial_spacetime_index_positions =
     354             :       detail::get_spatial_spacetime_index_positions<typename T2::index_list,
     355             :                                                     ArgsList2<Args2...>>();
     356             : 
     357             :   /// Whether or not the two operands have the same `TensorIndex`s in the same
     358             :   /// order (including concrete time indices)
     359           1 :   static constexpr bool ops_have_generic_indices_at_same_positions =
     360             :       generic_indices_at_same_positions<tmpl::list<Args1...>,
     361             :                                         tmpl::list<Args2...>>::value;
     362             : 
     363             :   // === Expression subtree properties ===
     364             :   /// The number of arithmetic tensor operations done in the subtree for the
     365             :   /// left operand
     366           1 :   static constexpr size_t num_ops_left_child = T1::num_ops_subtree;
     367             :   /// The number of arithmetic tensor operations done in the subtree for the
     368             :   /// right operand
     369           1 :   static constexpr size_t num_ops_right_child = T2::num_ops_subtree;
     370             :   // This helps ensure we are minimizing breadth in the overall tree when we
     371             :   // have addition (subtraction is not commutative)
     372             :   static_assert(Sign == -1 or num_ops_left_child >= num_ops_right_child,
     373             :                 "The left operand of an AddSub expression performing addition "
     374             :                 "should be a subtree with equal or more tensor operations than "
     375             :                 "the right operand's subtree.");
     376             :   /// The total number of arithmetic tensor operations done in this expression's
     377             :   /// whole subtree
     378           1 :   static constexpr size_t num_ops_subtree =
     379             :       num_ops_left_child + num_ops_right_child + 1;
     380             :   /// The height of this expression's node in the expression tree relative to
     381             :   /// the closest `TensorAsExpression` leaf in its subtree
     382           1 :   static constexpr size_t height_relative_to_closest_tensor_leaf_in_subtree =
     383             :       T2::height_relative_to_closest_tensor_leaf_in_subtree <=
     384             :               T1::height_relative_to_closest_tensor_leaf_in_subtree
     385             :           ? (T2::height_relative_to_closest_tensor_leaf_in_subtree !=
     386             :                      std::numeric_limits<size_t>::max()
     387             :                  ? T2::height_relative_to_closest_tensor_leaf_in_subtree + 1
     388             :                  : T2::height_relative_to_closest_tensor_leaf_in_subtree)
     389             :           : T1::height_relative_to_closest_tensor_leaf_in_subtree + 1;
     390             : 
     391             :   // === Properties for splitting up subexpressions along the primary path ===
     392             :   // These definitions only have meaning if this expression actually ends up
     393             :   // being along the primary path that is taken when evaluating the whole tree.
     394             :   // See documentation for `TensorExpression` for more details.
     395             :   /// If on the primary path, whether or not the expression is an ending point
     396             :   /// of a leg
     397           1 :   static constexpr bool is_primary_end = T1::is_primary_start;
     398             :   /// If on the primary path, this is the remaining number of arithmetic tensor
     399             :   /// operations that need to be done in the subtree of the child along the
     400             :   /// primary path, given that we will have already computed the whole subtree
     401             :   /// at the next lowest leg's starting point.
     402           1 :   static constexpr size_t num_ops_to_evaluate_primary_left_child =
     403             :       is_primary_end ? 0 : T1::num_ops_to_evaluate_primary_subtree;
     404             :   /// If on the primary path, this is the remaining number of arithmetic tensor
     405             :   /// operations that need to be done in the right operand's subtree. No
     406             :   /// splitting is currently done, so this is just `num_ops_right_child`.
     407           1 :   static constexpr size_t num_ops_to_evaluate_primary_right_child =
     408             :       num_ops_right_child;
     409             :   /// If on the primary path, this is the remaining number of arithmetic tensor
     410             :   /// operations that need to be done for this expression's subtree, given that
     411             :   /// we will have already computed the subtree at the next lowest leg's
     412             :   /// starting point
     413           1 :   static constexpr size_t num_ops_to_evaluate_primary_subtree =
     414             :       num_ops_to_evaluate_primary_left_child +
     415             :       num_ops_to_evaluate_primary_right_child + 1;
     416             :   /// If on the primary path, whether or not the expression is a starting point
     417             :   /// of a leg
     418           1 :   static constexpr bool is_primary_start =
     419             :       num_ops_to_evaluate_primary_subtree >=
     420             :       detail::max_num_ops_in_sub_expression<type>;
     421             :   /// When evaluating along a primary path, whether each operand's subtrees
     422             :   /// should be evaluated separately. Since `DataVector` expression runtime
     423             :   /// scales poorly with increased number of operations, evaluating the two
     424             :   /// expression subtrees separately like this is beneficial when at least one
     425             :   /// of the subtrees contains a large number of operations.
     426           1 :   static constexpr bool evaluate_children_separately =
     427             :       is_primary_start and (num_ops_to_evaluate_primary_left_child >=
     428             :                                 detail::max_num_ops_in_sub_expression<type> or
     429             :                             num_ops_to_evaluate_primary_right_child >=
     430             :                                 detail::max_num_ops_in_sub_expression<type>);
     431             :   /// If on the primary path, whether or not the expression's child along the
     432             :   /// primary path is a subtree that contains a starting point of a leg along
     433             :   /// the primary path
     434           1 :   static constexpr bool primary_child_subtree_contains_primary_start =
     435             :       T1::primary_subtree_contains_primary_start;
     436             :   /// If on the primary path, whether or not this subtree contains a starting
     437             :   /// point of a leg along the primary path
     438           1 :   static constexpr bool primary_subtree_contains_primary_start =
     439             :       is_primary_start or primary_child_subtree_contains_primary_start;
     440             : 
     441           0 :   AddSub(T1 t1, T2 t2) : t1_(std::move(t1)), t2_(std::move(t2)) {}
     442           0 :   ~AddSub() override = default;
     443             : 
     444             :   /// \brief Assert that the LHS tensor of the equation does not also appear in
     445             :   /// this expression's subtree
     446             :   template <typename LhsTensor>
     447           1 :   SPECTRE_ALWAYS_INLINE void assert_lhs_tensor_not_in_rhs_expression(
     448             :       const gsl::not_null<LhsTensor*> lhs_tensor) const {
     449             :     if constexpr (not std::is_base_of_v<MarkAsNumberAsExpression, T1>) {
     450             :       t1_.assert_lhs_tensor_not_in_rhs_expression(lhs_tensor);
     451             :     }
     452             :     if constexpr (not std::is_base_of_v<MarkAsNumberAsExpression, T2>) {
     453             :       t2_.assert_lhs_tensor_not_in_rhs_expression(lhs_tensor);
     454             :     }
     455             :   }
     456             : 
     457             :   /// \brief Assert that each instance of the LHS tensor in the RHS tensor
     458             :   /// expression uses the same generic index order that the LHS uses
     459             :   ///
     460             :   /// \tparam LhsTensorIndices the list of generic `TensorIndex`s of the LHS
     461             :   /// result `Tensor` being computed
     462             :   /// \param lhs_tensor the LHS result `Tensor` being computed
     463             :   template <typename LhsTensorIndices, typename LhsTensor>
     464           1 :   SPECTRE_ALWAYS_INLINE void assert_lhs_tensorindices_same_in_rhs(
     465             :       const gsl::not_null<LhsTensor*> lhs_tensor) const {
     466             :     if constexpr (not std::is_base_of_v<MarkAsNumberAsExpression, T1>) {
     467             :       t1_.template assert_lhs_tensorindices_same_in_rhs<LhsTensorIndices>(
     468             :           lhs_tensor);
     469             :     }
     470             :     if constexpr (not std::is_base_of_v<MarkAsNumberAsExpression, T2>) {
     471             :       t2_.template assert_lhs_tensorindices_same_in_rhs<LhsTensorIndices>(
     472             :           lhs_tensor);
     473             :     }
     474             :   }
     475             : 
     476             :   /// \brief Get the size of a component from a `Tensor` in this expression's
     477             :   /// subtree of the RHS `TensorExpression`
     478             :   ///
     479             :   /// \return the size of a component from a `Tensor` in this expression's
     480             :   /// subtree of the RHS `TensorExpression`
     481           1 :   SPECTRE_ALWAYS_INLINE size_t get_rhs_tensor_component_size() const {
     482             :     if constexpr (T1::height_relative_to_closest_tensor_leaf_in_subtree <=
     483             :                   T2::height_relative_to_closest_tensor_leaf_in_subtree) {
     484             :       return t1_.get_rhs_tensor_component_size();
     485             :     } else {
     486             :       return t2_.get_rhs_tensor_component_size();
     487             :     }
     488             :   }
     489             : 
     490             :   /// \brief Return the second operand's multi-index given the first operand's
     491             :   /// multi-index
     492             :   ///
     493             :   /// \param op1_multi_index the multi-index of the left operand
     494             :   /// \return the second operand's multi-index
     495             :   SPECTRE_ALWAYS_INLINE std::array<size_t, num_tensor_indices_op2>
     496           1 :   get_op2_multi_index(
     497             :       const std::array<size_t, num_tensor_indices>& op1_multi_index) const {
     498             :     if constexpr (ops_have_generic_indices_at_same_positions) {
     499             :       if constexpr (op1_spatial_spacetime_index_positions.size() != 0 or
     500             :                     op2_spatial_spacetime_index_positions.size() != 0) {
     501             :         // Operands have the same generic index order, but at least one of them
     502             :         // has at least one spacetime index where a spatial index has been used,
     503             :         // so we need to compute the 2nd operand's (possibly) shifted
     504             :         // multi-index values
     505             :         constexpr std::array<std::int32_t, num_tensor_indices>
     506             :             spatial_spacetime_index_transformation =
     507             :                 detail::spatial_spacetime_index_transformation_from_positions<
     508             :                     num_tensor_indices>(op1_spatial_spacetime_index_positions,
     509             :                                         op2_spatial_spacetime_index_positions);
     510             :         std::array<size_t, num_tensor_indices> op2_multi_index =
     511             :             op1_multi_index;
     512             :         for (size_t i = 0; i < num_tensor_indices; i++) {
     513             :           gsl::at(op2_multi_index, i) = static_cast<size_t>(
     514             :               static_cast<std::int32_t>(gsl::at(op2_multi_index, i)) +
     515             :               gsl::at(spatial_spacetime_index_transformation, i));
     516             :         }
     517             :         return op2_multi_index;
     518             :       } else {
     519             :         // Operands have the same generic index order and neither of them has
     520             :         // a spacetime index where a spatial index has been used, so
     521             :         // both operands have the same multi-index
     522             :         return op1_multi_index;
     523             :       }
     524             :     } else {
     525             :       if constexpr (op1_spatial_spacetime_index_positions.size() != 0 or
     526             :                     op2_spatial_spacetime_index_positions.size() != 0) {
     527             :         // Operands don't have the same generic index order and at least one of
     528             :         // them has at least one spacetime index where a spatial index has been
     529             :         // used, so we need to compute the 2nd operand's (possibly) shifted
     530             :         // multi-index values and reorder them with respect to the 2nd operand's
     531             :         // generic index order
     532             : 
     533             :         // The list of positions where generic spatial indices were used for
     534             :         // spacetime indices in the second operand, but rearranged in terms of
     535             :         // the first operand's generic index order.
     536             :         constexpr std::array<size_t,
     537             :                              op2_spatial_spacetime_index_positions.size()>
     538             :             transformed_op2_spatial_spacetime_index_positions = []() {
     539             :               std::array<size_t, op2_spatial_spacetime_index_positions.size()>
     540             :                   transformed_positions{};
     541             :               for (size_t i = 0;
     542             :                    i < op2_spatial_spacetime_index_positions.size(); i++) {
     543             :                 gsl::at(transformed_positions, i) =
     544             :                     gsl::at(operand_index_transformation,
     545             :                             gsl::at(op2_spatial_spacetime_index_positions, i));
     546             :               }
     547             :               return transformed_positions;
     548             :             }();
     549             : 
     550             :         // According to the transformed positions above, compute the value shift
     551             :         // needed to convert from multi-indices of the first operand to
     552             :         // multi-indices of the 2nd operand (with the generic index order of the
     553             :         // first)
     554             :         constexpr std::array<std::int32_t, num_tensor_indices>
     555             :             spatial_spacetime_index_transformation =
     556             :                 detail::spatial_spacetime_index_transformation_from_positions<
     557             :                     num_tensor_indices>(
     558             :                     op1_spatial_spacetime_index_positions,
     559             :                     transformed_op2_spatial_spacetime_index_positions);
     560             :         std::array<size_t, num_tensor_indices> op2_multi_index =
     561             :             op1_multi_index;
     562             :         for (size_t i = 0; i < num_tensor_indices; i++) {
     563             :           gsl::at(op2_multi_index, i) = static_cast<size_t>(
     564             :               static_cast<std::int32_t>(gsl::at(op2_multi_index, i)) +
     565             :               gsl::at(spatial_spacetime_index_transformation, i));
     566             :         }
     567             :         return transform_multi_index(op2_multi_index,
     568             :                                      operand_index_transformation);
     569             :       } else {
     570             :         // Operands don't have the same generic index order, but neither of them
     571             :         // has a spacetime index where a spatial index has been used, so we just
     572             :         // need to reorder the 2nd operand's multi_index according to its
     573             :         // generic index order
     574             :         return transform_multi_index(op1_multi_index,
     575             :                                      operand_index_transformation);
     576             :       }
     577             :     }
     578             :   }
     579             : 
     580             :   /// \brief Helper function for computing the sum of or difference between
     581             :   /// components at given multi-indices from both operands of the expression
     582             :   ///
     583             :   /// \details Both multi-index arguments must be ordered according to their
     584             :   /// operand's respective generic index ordering
     585             :   ///
     586             :   /// \param op1_multi_index the multi-index of the component of the first
     587             :   /// operand
     588             :   /// \param op2_multi_index the multi-index of the component of the second
     589             :   /// operand
     590             :   /// \return the sum of or difference between the two components' values
     591           1 :   SPECTRE_ALWAYS_INLINE decltype(auto) add_or_subtract(
     592             :       const std::array<size_t, num_tensor_indices>& op1_multi_index,
     593             :       const std::array<size_t, num_tensor_indices_op2>& op2_multi_index) const {
     594             :     if constexpr (Sign == 1) {
     595             :       return t1_.get(op1_multi_index) + t2_.get(op2_multi_index);
     596             :     } else {
     597             :       return t1_.get(op1_multi_index) - t2_.get(op2_multi_index);
     598             :     }
     599             :   }
     600             : 
     601             :   /// \brief Return the value of the component at the given multi-index of the
     602             :   /// tensor resulting from addition or subtraction
     603             :   ///
     604             :   /// \details One important detail to note about the type of the `AddSub`
     605             :   /// expression is that its two operands may have (i) different generic index
     606             :   /// orders, and/or (ii) different indices in their `index_list`s if where one
     607             :   /// operand uses a generic spatial index for a spacetime index, the other
     608             :   /// tensor may use that generic spatial index for a spatial index of the same
     609             :   /// dimension, valence, and frame. Therefore, there are four possible cases
     610             :   /// for an `AddSub` expression that are considered in the implementation:
     611             :   /// - same generic index order, spatial spacetime indices in expression
     612             :   /// - same generic index order, spatial spacetime indices not in expression
     613             :   /// - different generic index order, spatial spacetime indices in expression
     614             :   /// - different generic index order, spatial spacetime indices not in
     615             :   /// expression
     616             :   ///
     617             :   /// This means that for expressions where the generic index orders differ, a
     618             :   /// multi-index for a component of one operand is a (possible) rearrangement
     619             :   /// of the equivalent multi-index for a component in the other operand. This
     620             :   /// also means that for expressions where (at least once) a generic spatial
     621             :   /// index is used for a spacetime index, then, after accounting
     622             :   /// for potential reordering due to different generic index orders, a
     623             :   /// multi-index's values for a component of one operand are (possibly) shifted
     624             :   /// by one, compared to the multi-index's values for a component in the other
     625             :   /// operand.
     626             :   ///
     627             :   /// For example, given \f$R_{ij} + S_{ji}\f$, let \f$R\f$'s first index be
     628             :   /// a spacetime index, but \f$R\f$'s second index and both of \f$S\f$' indices
     629             :   /// be spatial indices. If \f$i = 2\f$ and \f$j = 0\f$, then when we compute
     630             :   /// \f$R_{20} + S_{02}\f$, the multi-index for \f$R_{20}\f$ is
     631             :   /// `{2 + 1, 0} = {3, 0}` (first value shifted because it is a spacetime
     632             :   /// index) and the multi-index for \f$S_{02}\f$ is `[0, 2]`. Because the first
     633             :   /// operand of an `AddSub` expresion propagates its generic index order and
     634             :   /// index list ( \ref SpacetimeIndex "TensorIndexType"s) as the `AddSub`'s own
     635             :   /// generic index order and index list, the `result_multi_index` is equivalent
     636             :   /// to the multi-index for the first operand. Thus, we need only compute the
     637             :   /// second operand's multi-index as a transformation of the first: reorder and
     638             :   /// shift the values of the first operand to compute the equivalent
     639             :   /// multi-index for the second operand.
     640             :   ///
     641             :   /// \param result_multi_index the multi-index of the component of the result
     642             :   /// tensor to retrieve
     643             :   /// \return the value of the component at `result_multi_index` in the result
     644             :   /// tensor
     645           1 :   SPECTRE_ALWAYS_INLINE decltype(auto) get(
     646             :       const std::array<size_t, num_tensor_indices>& result_multi_index) const {
     647             :     return add_or_subtract(result_multi_index,
     648             :                            get_op2_multi_index(result_multi_index));
     649             :   }
     650             : 
     651             :   /// \brief Helper for evaluating the LHS Tensor's result component at this
     652             :   /// subtree by evaluating the two operand's subtrees separately and adding or
     653             :   /// subtracting them
     654             :   ///
     655             :   /// \details
     656             :   /// The left and right operands' subtrees are evaluated successively with
     657             :   /// two separate assignments to the LHS result component. Since `DataVector`
     658             :   /// expression runtime scales poorly with increased number of operations,
     659             :   /// evaluating the two expression subtrees separately like this is beneficial
     660             :   /// when at least one of the subtrees contains a large number of operations.
     661             :   /// Instead of evaluating a larger expression with their combined total number
     662             :   /// of operations, we evaluate two smaller ones.
     663             :   ///
     664             :   /// This function also differs from `add_or_subtract` in that it takes into
     665             :   /// account whether we have already computed part of the result component at a
     666             :   /// lower subtree. In recursively computing this sum/difference, the current
     667             :   /// result component will be substituted in for the most recent (highest)
     668             :   /// subtree below it that has already been evaluated.
     669             :   ///
     670             :   /// \param result_component the LHS tensor component to evaluate
     671             :   /// \param op1_multi_index the multi-index of the component of the first
     672             :   /// operand of the sum or difference to evaluate
     673             :   /// \param op2_multi_index the multi-index of the component of the second
     674             :   /// operand of the sum or difference to evaluate
     675             :   template <typename ResultType>
     676           1 :   SPECTRE_ALWAYS_INLINE void add_or_subtract_primary_children(
     677             :       ResultType& result_component,
     678             :       const std::array<size_t, num_tensor_indices>& op1_multi_index,
     679             :       const std::array<size_t, num_tensor_indices_op2>& op2_multi_index) const {
     680             :     if constexpr (Sign == 1) {
     681             :       // We're performing addition
     682             :       if constexpr (is_primary_end) {
     683             :         (void)op1_multi_index;
     684             :         // We've already computed the whole child subtree on the primary path,
     685             :         // so just add the result of the other child's subtree to the current
     686             :         // result
     687             :         result_component += t2_.get(op2_multi_index);
     688             :       } else {
     689             :         // We haven't yet evaluated the whole subtree of the primary child, so
     690             :         // first assign the result component to be the result of computing the
     691             :         // primary child's subtree
     692             :         result_component =
     693             :             t1_.template get_primary(result_component, op1_multi_index);
     694             :         // Now that the primary child's subtree has been computed, add the
     695             :         // result of evaluating the other child's subtree to the current result
     696             :         result_component += t2_.get(op2_multi_index);
     697             :       }
     698             :     } else {
     699             :       // We're performing subtraction
     700             :       if constexpr (is_primary_end) {
     701             :         (void)op1_multi_index;
     702             :         // We've already computed the whole child subtree on the primary path,
     703             :         // so just subtract the result of the other child's subtree from the
     704             :         // current result
     705             :         result_component -= t2_.get(op2_multi_index);
     706             :       } else {
     707             :         // We haven't yet evaluated the whole subtree of the primary child, so
     708             :         // first assign the result component to be the result of computing the
     709             :         // primary child's subtree
     710             :         result_component =
     711             :             t1_.template get_primary(result_component, op1_multi_index);
     712             :         // Now that the primary child's subtree has been computed, subtract the
     713             :         // result of evaluating the other child's subtree from the current
     714             :         // result
     715             :         result_component -= t2_.get(op2_multi_index);
     716             :       }
     717             :     }
     718             :   }
     719             : 
     720             :   /// \brief Evaluate the LHS Tensor's result component at this subtree by
     721             :   /// evaluating the two operand's subtrees separately and adding or subtracting
     722             :   /// them
     723             :   ///
     724             :   /// \details
     725             :   /// See `add_or_subtract_primary_children` for more details
     726             :   ///
     727             :   /// \param result_component the LHS tensor component to evaluate
     728             :   /// \param result_multi_index the multi-index of the component of the result
     729             :   /// tensor to evaluate
     730             :   template <typename ResultType>
     731           1 :   SPECTRE_ALWAYS_INLINE void evaluate_primary_children(
     732             :       ResultType& result_component,
     733             :       const std::array<size_t, num_tensor_indices>& result_multi_index) const {
     734             :     add_or_subtract_primary_children(result_component, result_multi_index,
     735             :                                      get_op2_multi_index(result_multi_index));
     736             :   }
     737             : 
     738             :   /// \brief Helper function for returning the sum of or difference between
     739             :   /// components at given multi-indices from both operands of the expression
     740             :   ///
     741             :   /// \details
     742             :   /// This function differs from `add_or_subtract` in that it takes into account
     743             :   /// whether we have already computed part of the result component at a lower
     744             :   /// subtree. In recursively computing this sum/difference, the current result
     745             :   /// component will be substituted in for the most recent (highest) subtree
     746             :   /// below it that has already been evaluated.
     747             :   ///
     748             :   /// \param result_component the LHS tensor component to evaluate
     749             :   /// \param op1_multi_index the multi-index of the component of the first
     750             :   /// operand
     751             :   /// \param op2_multi_index the multi-index of the component of the second
     752             :   /// operand
     753             :   /// \return the sum of or difference between the two components' values
     754             :   template <typename ResultType>
     755           1 :   SPECTRE_ALWAYS_INLINE decltype(auto) add_or_subtract_primary(
     756             :       const ResultType& result_component,
     757             :       const std::array<size_t, num_tensor_indices>& op1_multi_index,
     758             :       const std::array<size_t, num_tensor_indices_op2>& op2_multi_index) const {
     759             :     if constexpr (Sign == 1) {
     760             :       // We're performing addition
     761             :       if constexpr (is_primary_end) {
     762             :         (void)op1_multi_index;
     763             :         // We've already computed the whole child subtree on the primary path,
     764             :         // so just add the result of the other child's subtree to the current
     765             :         // result
     766             :         return result_component + t2_.get(op2_multi_index);
     767             :       } else {
     768             :         // We haven't yet evaluated the whole subtree for this expression, so
     769             :         // return the sum of the results of the two operands' subtrees
     770             :         return t1_.template get_primary(result_component, op1_multi_index) +
     771             :                t2_.get(op2_multi_index);
     772             :       }
     773             :     } else {
     774             :       // We're performing subtraction
     775             :       if constexpr (is_primary_end) {
     776             :         (void)op1_multi_index;
     777             :         // We've already computed the whole child subtree on the primary path,
     778             :         // so just subtract the result of the other child's subtree from the
     779             :         // current result
     780             :         return result_component - t2_.get(op2_multi_index);
     781             :       } else {
     782             :         // We haven't yet evaluated the whole subtree for this expression, so
     783             :         // return the difference between the results of the two operands'
     784             :         // subtrees
     785             :         return t1_.template get_primary(result_component, op1_multi_index) -
     786             :                t2_.get(op2_multi_index);
     787             :       }
     788             :     }
     789             :   }
     790             : 
     791             :   /// \brief Return the value of the component at the given multi-index of the
     792             :   /// tensor resulting from addition or subtraction
     793             :   ///
     794             :   /// \details
     795             :   /// This function differs from `get` in that it takes into account whether we
     796             :   /// have already computed part of the result component at a lower subtree.
     797             :   /// In recursively computing this sum/difference, the current result component
     798             :   /// will be substituted in for the most recent (highest) subtree below it that
     799             :   /// has already been evaluated.
     800             :   ///
     801             :   /// \param result_component the LHS tensor component to evaluate
     802             :   /// \param result_multi_index the multi-index of the component of the result
     803             :   /// tensor to retrieve
     804             :   /// \return the value of the component at `result_multi_index` in the result
     805             :   /// tensor
     806             :   template <typename ResultType>
     807           1 :   SPECTRE_ALWAYS_INLINE decltype(auto) get_primary(
     808             :       const ResultType& result_component,
     809             :       const std::array<size_t, num_tensor_indices>& result_multi_index) const {
     810             :     return add_or_subtract_primary(result_component, result_multi_index,
     811             :                                    get_op2_multi_index(result_multi_index));
     812             :   }
     813             : 
     814             :   /// \brief Successively evaluate the LHS Tensor's result component at each
     815             :   /// leg in this expression's subtree
     816             :   ///
     817             :   /// \details
     818             :   /// This function takes into account whether we have already computed part of
     819             :   /// the result component at a lower subtree. In recursively computing this
     820             :   /// sum/difference, the current result component will be substituted in for
     821             :   /// the most recent (highest) subtree below it that has already been
     822             :   /// evaluated.
     823             :   ///
     824             :   /// \param result_component the LHS tensor component to evaluate
     825             :   /// \param result_multi_index the multi-index of the component of the result
     826             :   /// tensor to evaluate
     827             :   template <typename ResultType>
     828           1 :   SPECTRE_ALWAYS_INLINE void evaluate_primary_subtree(
     829             :       ResultType& result_component,
     830             :       const std::array<size_t, num_tensor_indices>& result_multi_index) const {
     831             :     if constexpr (primary_child_subtree_contains_primary_start) {
     832             :       // The primary child's subtree contains at least one leg, so recurse down
     833             :       // and evaluate that first
     834             :       t1_.template evaluate_primary_subtree(result_component,
     835             :                                             result_multi_index);
     836             :     }
     837             : 
     838             :     if constexpr (is_primary_start) {
     839             :       // We want to evaluate the subtree for this expression
     840             :       if constexpr (evaluate_children_separately) {
     841             :         // Evaluate operand's subtrees separately
     842             :         evaluate_primary_children(result_component, result_multi_index);
     843             :       } else {
     844             :         // Evaluate whole subtree as one expression
     845             :         result_component = get_primary(result_component, result_multi_index);
     846             :       }
     847             :     }
     848             :   }
     849             : 
     850             :  private:
     851             :   /// Left operand expression
     852           1 :   T1 t1_;
     853             :   /// Right operand expression
     854           1 :   T2 t2_;
     855             : };
     856             : }  // namespace tenex
     857             : 
     858             : /*!
     859             :  * \ingroup TensorExpressionsGroup
     860             :  */
     861             : template <typename T1, typename T2, typename X1, typename X2, typename Symm1,
     862             :           typename Symm2, typename IndexList1, typename IndexList2,
     863             :           typename Args1, typename Args2>
     864           0 : SPECTRE_ALWAYS_INLINE auto operator+(
     865             :     const TensorExpression<T1, X1, Symm1, IndexList1, Args1>& t1,
     866             :     const TensorExpression<T2, X2, Symm2, IndexList2, Args2>& t2) {
     867             :   using op1_generic_indices =
     868             :       typename tenex::detail::remove_time_indices<Args1>::type;
     869             :   using op2_generic_indices =
     870             :       typename tenex::detail::remove_time_indices<Args2>::type;
     871             :   static_assert(tmpl::size<op1_generic_indices>::value ==
     872             :                     tmpl::size<op2_generic_indices>::value,
     873             :                 "Tensor addition is only possible when the same number of "
     874             :                 "generic indices are used with both operands.");
     875             :   static_assert(
     876             :       tmpl::equal_members<op1_generic_indices, op2_generic_indices>::value,
     877             :       "The generic indices when adding two tensors must be equal. This error "
     878             :       "occurs from expressions like R(ti::a, ti::b) + S(ti::c, ti::a)");
     879             :   if constexpr (T1::num_ops_subtree >= T2::num_ops_subtree) {
     880             :     return tenex::AddSub<T1, T2, Args1, Args2, 1>(~t1, ~t2);
     881             :   } else {
     882             :     return tenex::AddSub<T2, T1, Args2, Args1, 1>(~t2, ~t1);
     883             :   }
     884             : }
     885             : 
     886             : /// @{
     887             : /// \ingroup TensorExpressionsGroup
     888             : /// \brief Returns the tensor expression representing the sum of a tensor
     889             : /// expression and a number
     890             : ///
     891             : /// \details
     892             : /// The tensor expression operand must represent an expression that, when
     893             : /// evaluated, would be a rank 0 tensor. For example, if `R` and `S` are
     894             : /// Tensors, here is a non-exhaustive list of some of the acceptable forms that
     895             : /// the tensor expression operand could take:
     896             : /// - `R()`
     897             : /// - `R(ti::A, ti::a)`
     898             : /// - `(R(ti::A, ti::B) * S(ti::a, ti::b))`
     899             : /// - `R(ti::t, ti::t)`
     900             : ///
     901             : /// \param t the tensor expression operand of the sum
     902             : /// \param number the numeric operand of the sum
     903             : /// \return the tensor expression representing the sum of the tensor expression
     904             : /// and the number
     905             : template <typename T, typename X, typename Symm, typename IndexList,
     906             :           typename... Args, typename N,
     907             :           Requires<std::is_arithmetic_v<N>> = nullptr>
     908           1 : SPECTRE_ALWAYS_INLINE auto operator+(
     909             :     const TensorExpression<T, X, Symm, IndexList, tmpl::list<Args...>>& t,
     910             :     const N number) {
     911             :   static_assert(
     912             :       (... and tt::is_time_index<Args>::value),
     913             :       "Can only add a number to a tensor expression that evaluates to a rank 0"
     914             :       "tensor.");
     915             :   return t + tenex::NumberAsExpression(number);
     916             : }
     917             : template <typename T, typename X, typename Symm, typename IndexList,
     918             :           typename... Args, typename N,
     919             :           Requires<std::is_arithmetic_v<N>> = nullptr>
     920           1 : SPECTRE_ALWAYS_INLINE auto operator+(
     921             :     const N number,
     922             :     const TensorExpression<T, X, Symm, IndexList, tmpl::list<Args...>>& t) {
     923             :   static_assert(
     924             :       (... and tt::is_time_index<Args>::value),
     925             :       "Can only add a number to a tensor expression that evaluates to a rank 0"
     926             :       "tensor.");
     927             :   return t + tenex::NumberAsExpression(number);
     928             : }
     929             : template <typename T, typename X, typename Symm, typename IndexList,
     930             :           typename... Args, typename N>
     931           1 : SPECTRE_ALWAYS_INLINE auto operator+(
     932             :     const TensorExpression<T, X, Symm, IndexList, tmpl::list<Args...>>& t,
     933             :     const std::complex<N>& number) {
     934             :   static_assert(
     935             :       (... and tt::is_time_index<Args>::value),
     936             :       "Can only add a number to a tensor expression that evaluates to a rank 0"
     937             :       "tensor.");
     938             :   return t + tenex::NumberAsExpression(number);
     939             : }
     940             : template <typename T, typename X, typename Symm, typename IndexList,
     941             :           typename... Args, typename N>
     942           1 : SPECTRE_ALWAYS_INLINE auto operator+(
     943             :     const std::complex<N>& number,
     944             :     const TensorExpression<T, X, Symm, IndexList, tmpl::list<Args...>>& t) {
     945             :   static_assert(
     946             :       (... and tt::is_time_index<Args>::value),
     947             :       "Can only add a number to a tensor expression that evaluates to a rank 0"
     948             :       "tensor.");
     949             :   return t + tenex::NumberAsExpression(number);
     950             : }
     951             : /// @}
     952             : 
     953             : /*!
     954             :  * \ingroup TensorExpressionsGroup
     955             :  */
     956             : template <typename T1, typename T2, typename X1, typename X2, typename Symm1,
     957             :           typename Symm2, typename IndexList1, typename IndexList2,
     958             :           typename Args1, typename Args2>
     959           0 : SPECTRE_ALWAYS_INLINE auto operator-(
     960             :     const TensorExpression<T1, X1, Symm1, IndexList1, Args1>& t1,
     961             :     const TensorExpression<T2, X2, Symm2, IndexList2, Args2>& t2) {
     962             :   using op1_generic_indices =
     963             :       typename tenex::detail::remove_time_indices<Args1>::type;
     964             :   using op2_generic_indices =
     965             :       typename tenex::detail::remove_time_indices<Args2>::type;
     966             :   static_assert(tmpl::size<op1_generic_indices>::value ==
     967             :                     tmpl::size<op2_generic_indices>::value,
     968             :                 "Tensor subtraction is only possible when the same number of "
     969             :                 "generic indices are used with both operands.");
     970             :   static_assert(
     971             :       tmpl::equal_members<op1_generic_indices, op2_generic_indices>::value,
     972             :       "The generic indices when subtracting two tensors must be equal. This "
     973             :       "error occurs from expressions like R(ti::a, ti::b) - S(ti::c, ti::a)");
     974             :   return tenex::AddSub<T1, T2, Args1, Args2, -1>(~t1, ~t2);
     975             : }
     976             : 
     977             : /// @{
     978             : /// \ingroup TensorExpressionsGroup
     979             : /// \brief Returns the tensor expression representing the difference of a tensor
     980             : /// expression and a number
     981             : ///
     982             : /// \details
     983             : /// The tensor expression operand must represent an expression that, when
     984             : /// evaluated, would be a rank 0 tensor. For example, if `R` and `S` are
     985             : /// Tensors, here is a non-exhaustive list of some of the acceptable forms that
     986             : /// the tensor expression operand could take:
     987             : /// - `R()`
     988             : /// - `R(ti::A, ti::a)`
     989             : /// - `(R(ti::A, ti::B) * S(ti::a, ti::b))`
     990             : /// - `R(ti::t, ti::t)`
     991             : ///
     992             : /// \param t the tensor expression operand of the difference
     993             : /// \param number the numeric operand of the difference
     994             : /// \return the tensor expression representing the difference of the tensor
     995             : /// expression and the number
     996             : template <typename T, typename X, typename Symm, typename IndexList,
     997             :           typename... Args, typename N,
     998             :           Requires<std::is_arithmetic_v<N>> = nullptr>
     999           1 : SPECTRE_ALWAYS_INLINE auto operator-(
    1000             :     const TensorExpression<T, X, Symm, IndexList, tmpl::list<Args...>>& t,
    1001             :     const N number) {
    1002             :   static_assert(
    1003             :       (... and tt::is_time_index<Args>::value),
    1004             :       "Can only subtract a number from a tensor expression that evaluates to a "
    1005             :       "rank 0 tensor.");
    1006             :   return t + tenex::NumberAsExpression(-number);
    1007             : }
    1008             : template <typename T, typename X, typename Symm, typename IndexList,
    1009             :           typename... Args, typename N,
    1010             :           Requires<std::is_arithmetic_v<N>> = nullptr>
    1011           1 : SPECTRE_ALWAYS_INLINE auto operator-(
    1012             :     const N number,
    1013             :     const TensorExpression<T, X, Symm, IndexList, tmpl::list<Args...>>& t) {
    1014             :   static_assert(
    1015             :       (... and tt::is_time_index<Args>::value),
    1016             :       "Can only subtract a number from a tensor expression that evaluates to a "
    1017             :       "rank 0 tensor.");
    1018             :   return tenex::NumberAsExpression(number) - t;
    1019             : }
    1020             : template <typename T, typename X, typename Symm, typename IndexList,
    1021             :           typename... Args, typename N>
    1022           1 : SPECTRE_ALWAYS_INLINE auto operator-(
    1023             :     const TensorExpression<T, X, Symm, IndexList, tmpl::list<Args...>>& t,
    1024             :     const std::complex<N>& number) {
    1025             :   static_assert(
    1026             :       (... and tt::is_time_index<Args>::value),
    1027             :       "Can only subtract a number from a tensor expression that evaluates to a "
    1028             :       "rank 0 tensor.");
    1029             :   return t + tenex::NumberAsExpression(-number);
    1030             : }
    1031             : template <typename T, typename X, typename Symm, typename IndexList,
    1032             :           typename... Args, typename N>
    1033           1 : SPECTRE_ALWAYS_INLINE auto operator-(
    1034             :     const std::complex<N>& number,
    1035             :     const TensorExpression<T, X, Symm, IndexList, tmpl::list<Args...>>& t) {
    1036             :   static_assert(
    1037             :       (... and tt::is_time_index<Args>::value),
    1038             :       "Can only subtract a number from a tensor expression that evaluates to a "
    1039             :       "rank 0 tensor.");
    1040             :   return tenex::NumberAsExpression(number) - t;
    1041             : }
    1042             : /// @}

Generated by: LCOV version 1.14