Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <cstddef>
       8             : #include <tuple>
       9             : #include <utility>
      10             : 
      11             : #include "DataStructures/DataBox/PrefixHelpers.hpp"
      12             : #include "DataStructures/Tags/TempTensor.hpp"
      13             : #include "DataStructures/Tensor/IndexType.hpp"
      14             : #include "DataStructures/Tensor/Tensor.hpp"
      15             : #include "DataStructures/Tensor/TypeAliases.hpp"
      16             : #include "DataStructures/Variables.hpp"
      17             : #include "Utilities/Gsl.hpp"
      18             : #include "Utilities/Requires.hpp"
      19             : #include "Utilities/TMPL.hpp"
      20             : #include "Utilities/TypeTraits/IsA.hpp"
      21             : 
      22             : namespace Spectral {
      23             : /*!
      24             :  * \brief A compile-time function to determine parity, with respect to the
      25             :  * \f$x\f$ coordinate in an axisymmetric spacetime, of tensors in a Variables
      26             :  *
      27             :  * \details In a \f$d\f$-dimensional space where fields are regular everywhere
      28             :  * (particularly at the axis), we utilize the fact that smoothness corresponds
      29             :  * to a well-defined Taylor expansion in all \f$d\f$ coordinates. Viewing the
      30             :  * components of a rank \f$(p, q)\f$ tensor as smooth functions of \f$(x, y,
      31             :  * \ldots, z)\f$, consider a reflection transformation \f$x \to -x\f$. As a
      32             :  * coordinate transformation, the Jacobian will look like \f$\mathrm{diag}(-1,
      33             :  * 1, 1, \ldots)\f$.
      34             :  *
      35             :  * If any of the \f$(p, q)\f$ indices correspond to the \f$x\f$-coordinate for
      36             :  * this particular component, there is a \f$-1\f$ factor associated with each
      37             :  * \f$x\f$-index, i.e., the Jacobian will look like \f$(-1)^{n_x}\f$, where
      38             :  * \f$n_x\f$ is the number of \f$x\f$-indices of the component.
      39             :  *
      40             :  * If this space is symmetric under this reflection, as it is for axisymmetry
      41             :  * about the \f$y\f$-axis in 3 dimensions, components can only represent
      42             :  * functions invariant under this transformation.
      43             :  *
      44             :  * This means for an even number of \f$x\f$-indices, the Taylor expansion must
      45             :  * only have even powers of \f$x\f$, and for an odd number of \f$x\f$-indices,
      46             :  * the expansion must have purely odd powers of \f$x\f$.
      47             :  *
      48             :  * The primary use-case for this is axisymmetric Cartoon simulations, where
      49             :  * elements touching the symmetry axis require ZernikeB1 bases for numerical
      50             :  * stability. These bases require knowledge of component parity for certain
      51             :  * operations (differentiation, interpolation, etc.).
      52             :  *
      53             :  * The returned array, alternating values for even/odd, stores the number of
      54             :  * next components with the same parity. When there are any neighboring parities
      55             :  * (i.e. not the worst-case scenario) the array will be padded with zeros to
      56             :  * be the correct size. The first index may be 0 if the first component is
      57             :  * odd, but any following zeros are guaranteed to only have zeros following it.
      58             :  * The first returned `size_t` is the number of even components, while the
      59             :  * second is the number of odd components.
      60             :  */
      61             : template <typename VariablesTags>
      62             : constexpr std::tuple<
      63             :     std::array<size_t,
      64             :                Variables<VariablesTags>::number_of_independent_components + 1>,
      65             :     size_t, size_t>
      66           1 : compute_parity_list() {
      67             :   const size_t N =
      68             :       Variables<VariablesTags>::number_of_independent_components + 1;
      69             :   std::array<size_t, N> parity_run_lengths{};
      70             : 
      71             :   const auto is_x_coordinate_index = [](const IndexType index_type,
      72             :                                         const size_t index_value) {
      73             :     return (index_type == IndexType::Spacetime and index_value == 1) or
      74             :            (index_type == IndexType::Spatial and index_value == 0);
      75             :   };
      76             : 
      77             :   size_t run_index = 0;
      78             :   bool current_parity_is_even = true;
      79             :   tmpl::for_each<VariablesTags>([&parity_run_lengths, &run_index,
      80             :                                  &current_parity_is_even,
      81             :                                  &is_x_coordinate_index]<typename TensorTag>(
      82             :                                     tmpl::type_<TensorTag> /*meta*/) {
      83             :     using tensor_type = typename TensorTag::type;
      84             :     constexpr auto index_types = tensor_type::index_types();
      85             :     constexpr size_t tensor_size = tensor_type::size();
      86             : 
      87             :     for (size_t component_index = 0; component_index < tensor_size;
      88             :          ++component_index) {
      89             :       const auto tensor_index = tensor_type::get_tensor_index(component_index);
      90             :       size_t x_coordinate_count = 0;
      91             :       for (size_t index_position = 0; index_position < index_types.size();
      92             :            ++index_position) {
      93             :         if (is_x_coordinate_index(gsl::at(index_types, index_position),
      94             :                                   gsl::at(tensor_index, index_position))) {
      95             :           ++x_coordinate_count;
      96             :         }
      97             :       }
      98             :       // If current parity doesn't match last
      99             :       const bool component_is_even = (x_coordinate_count % 2 == 0);
     100             :       if (component_is_even != current_parity_is_even) {
     101             :         ++run_index;
     102             :         current_parity_is_even = !current_parity_is_even;
     103             :       }
     104             :       gsl::at(parity_run_lengths, run_index) += 1;
     105             :     }
     106             :   });
     107             : 
     108             :   const auto [num_even, num_odd] = [&parity_run_lengths]<size_t... Is>(
     109             :                                        std::index_sequence<Is...> /*meta*/) {
     110             :     std::size_t even_count = ((Is % 2 == 0 ? parity_run_lengths[Is] : 0) + ...);
     111             :     std::size_t odd_count = ((Is % 2 == 1 ? parity_run_lengths[Is] : 0) + ...);
     112             :     return std::pair{even_count, odd_count};
     113             :   }(std::make_index_sequence<N>{});
     114             : 
     115             :   return {parity_run_lengths, num_even, num_odd};
     116             : }
     117             : 
     118             : /*!
     119             :  * \brief A compile-time function to determine parity, with respect to the
     120             :  * \f$x\f$ coordinate in an axisymmetric spacetime, of a tensor.
     121             :  *
     122             :  * \see `compute_parity_list(Variables)`
     123             :  */
     124             : template <typename TensorType>
     125             :   requires(tt::is_a_v<Tensor, TensorType>)
     126             : constexpr std::tuple<std::array<size_t, TensorType::size() + 1>, size_t, size_t>
     127           1 : compute_parity_list() {
     128             :   using tensor_type = Tensor<DataVector, typename TensorType::symmetry,
     129             :                              typename TensorType::index_list>;
     130             :   using vars_list = ::Tags::convert_to_temp_tensors<tmpl::list<tensor_type>, 0>;
     131             :   return compute_parity_list<vars_list>();
     132             : }
     133             : }  // namespace Spectral