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 "NumericalAlgorithms/Spectral/Parity.hpp"
18 : #include "Utilities/Gsl.hpp"
19 : #include "Utilities/Requires.hpp"
20 : #include "Utilities/TMPL.hpp"
21 : #include "Utilities/TypeTraits/IsA.hpp"
22 :
23 : namespace Spectral {
24 : /*!
25 : * \brief A compile-time function to determine parity, with respect to the
26 : * \f$x\f$ coordinate in an axisymmetric spacetime, of tensors in a Variables
27 : *
28 : * \details In a \f$d\f$-dimensional space where fields are regular everywhere
29 : * (particularly at the axis), we utilize the fact that smoothness corresponds
30 : * to a well-defined Taylor expansion in all \f$d\f$ coordinates. Viewing the
31 : * components of a rank \f$(p, q)\f$ tensor as smooth functions of \f$(x, y,
32 : * \ldots, z)\f$, consider a reflection transformation \f$x \to -x\f$. As a
33 : * coordinate transformation, the Jacobian will look like \f$\mathrm{diag}(-1,
34 : * 1, 1, \ldots)\f$.
35 : *
36 : * If any of the \f$(p, q)\f$ indices correspond to the \f$x\f$-coordinate for
37 : * this particular component, there is a \f$-1\f$ factor associated with each
38 : * \f$x\f$-index, i.e., the Jacobian will look like \f$(-1)^{n_x}\f$, where
39 : * \f$n_x\f$ is the number of \f$x\f$-indices of the component.
40 : *
41 : * If this space is symmetric under this reflection, as it is for axisymmetry
42 : * about the \f$y\f$-axis in 3 dimensions, components can only represent
43 : * functions invariant under this transformation.
44 : *
45 : * This means for an even number of \f$x\f$-indices, the Taylor expansion must
46 : * only have even powers of \f$x\f$, and for an odd number of \f$x\f$-indices,
47 : * the expansion must have purely odd powers of \f$x\f$.
48 : *
49 : * The primary use-case for this is axisymmetric Cartoon simulations, where
50 : * elements touching the symmetry axis require ZernikeB1 bases for numerical
51 : * stability. These bases require knowledge of component parity for certain
52 : * operations (differentiation, interpolation, etc.).
53 : *
54 : * The returned array, alternating values for even/odd, stores the number of
55 : * next components with the same parity. When there are any neighboring parities
56 : * (i.e. not the worst-case scenario) the array will be padded with zeros to
57 : * be the correct size. The first index may be 0 if the first component is
58 : * odd, but any following zeros are guaranteed to only have zeros following it.
59 : * The first returned `size_t` is the number of even components, while the
60 : * second is the number of odd components.
61 : */
62 : template <typename VariablesTags>
63 : constexpr std::tuple<
64 : std::array<size_t,
65 : Variables<VariablesTags>::number_of_independent_components + 1>,
66 : size_t, size_t>
67 1 : compute_parity_list() {
68 : const size_t N =
69 : Variables<VariablesTags>::number_of_independent_components + 1;
70 : std::array<size_t, N> parity_run_lengths{};
71 :
72 : const auto is_x_coordinate_index = [](const IndexType index_type,
73 : const size_t index_value) {
74 : return (index_type == IndexType::Spacetime and index_value == 1) or
75 : (index_type == IndexType::Spatial and index_value == 0);
76 : };
77 :
78 : size_t run_index = 0;
79 : bool current_parity_is_even = true;
80 : tmpl::for_each<VariablesTags>([&parity_run_lengths, &run_index,
81 : ¤t_parity_is_even,
82 : &is_x_coordinate_index]<typename TensorTag>(
83 : tmpl::type_<TensorTag> /*meta*/) {
84 : using tensor_type = typename TensorTag::type;
85 : constexpr auto index_types = tensor_type::index_types();
86 : constexpr size_t tensor_size = tensor_type::size();
87 :
88 : for (size_t component_index = 0; component_index < tensor_size;
89 : ++component_index) {
90 : const auto tensor_index = tensor_type::get_tensor_index(component_index);
91 : size_t x_coordinate_count = 0;
92 : for (size_t index_position = 0; index_position < index_types.size();
93 : ++index_position) {
94 : if (is_x_coordinate_index(gsl::at(index_types, index_position),
95 : gsl::at(tensor_index, index_position))) {
96 : ++x_coordinate_count;
97 : }
98 : }
99 : // If current parity doesn't match last
100 : const bool component_is_even = (x_coordinate_count % 2 == 0);
101 : if (component_is_even != current_parity_is_even) {
102 : ++run_index;
103 : current_parity_is_even = !current_parity_is_even;
104 : }
105 : gsl::at(parity_run_lengths, run_index) += 1;
106 : }
107 : });
108 :
109 : const auto [num_even, num_odd] = [&parity_run_lengths]<size_t... Is>(
110 : std::index_sequence<Is...> /*meta*/) {
111 : std::size_t even_count = ((Is % 2 == 0 ? parity_run_lengths[Is] : 0) + ...);
112 : std::size_t odd_count = ((Is % 2 == 1 ? parity_run_lengths[Is] : 0) + ...);
113 : return std::pair{even_count, odd_count};
114 : }(std::make_index_sequence<N>{});
115 :
116 : return {parity_run_lengths, num_even, num_odd};
117 : }
118 :
119 : /*!
120 : * \brief A compile-time function to determine parity, with respect to the
121 : * \f$x\f$ coordinate in an axisymmetric spacetime, of a tensor.
122 : *
123 : * \see `compute_parity_list(Variables)`
124 : */
125 : template <typename TensorType>
126 : requires(tt::is_a_v<Tensor, TensorType>)
127 : constexpr std::tuple<std::array<size_t, TensorType::size() + 1>, size_t, size_t>
128 1 : compute_parity_list() {
129 : using tensor_type = Tensor<DataVector, typename TensorType::symmetry,
130 : typename TensorType::index_list>;
131 : using vars_list = ::Tags::convert_to_temp_tensors<tmpl::list<tensor_type>, 0>;
132 : return compute_parity_list<vars_list>();
133 : }
134 : /*!
135 : * \brief Returns a compile-time array mapping each component index of
136 : * `TensorType` to its `Parity` (Even or Odd) based on the number of
137 : * \f$x\f$-coordinate indices, for use in an axisymmetric spacetime with the
138 : * Cartoon method.
139 : *
140 : * A component is `Parity::Even` if it has an even number of
141 : * \f$x\f$-coordinate indices, or `Parity::Odd` otherwise.
142 : *
143 : * \see `compute_parity_list`
144 : */
145 : template <typename TensorType>
146 : requires(tt::is_a_v<Tensor, TensorType>)
147 1 : constexpr std::array<Parity, TensorType::size()> make_component_parity_array() {
148 : constexpr auto parity_info = compute_parity_list<TensorType>();
149 : constexpr auto parity_list = std::get<0>(parity_info);
150 : constexpr size_t N = TensorType::size();
151 : std::array<Parity, N> result{};
152 : size_t component = 0;
153 : bool is_even = true;
154 : for (size_t i = 0; component < N; ++i) {
155 : const size_t seg_size = parity_list[i];
156 : if (seg_size == 0) {
157 : if (is_even) {
158 : is_even = false;
159 : continue;
160 : } else {
161 : break;
162 : }
163 : }
164 : for (size_t k = 0; k < seg_size; ++k, ++component) {
165 : result[component] = is_even ? Parity::Even : Parity::Odd;
166 : }
167 : is_even = not is_even;
168 : }
169 : return result;
170 : }
171 : } // namespace Spectral
|
