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 of tensors in a
26 : * Variables in an axisymmetric spacetime.
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 above argument holds for operations related to the `ZernikeB1` basis,
50 : * which is used only in the \f$x\f$-direction of cartoon
51 : * simulations.
52 : *
53 : * The primary use-case is axisymmetric Cartoon simulations, where
54 : * elements touching the symmetry axis require ZernikeB1 in the \f$x\f$
55 : * direction bases for numerical stability. These bases require knowledge of
56 : * component parity for certain operations (differentiation, interpolation,
57 : * etc.). However, the true axial symmetry maps \f$(x, y, z) \to (-x, y,
58 : * -z)\f$. Both the \f$x\f$-direction and the \f$z\f$-direction flip sign,
59 : * contributing \f$(-1)^{n_x + n_z}\f$ per component. This is relevant for
60 : * operations with the `HalfFourier` basis when used in an axisymmetric
61 : * spacetime. Use `IncludeZ = true` for HalfFourier-based differentiation in
62 : * the azimuthal direction.
63 : *
64 : * The returned array, alternating values for even/odd, stores the number of
65 : * next components with the same parity. When there are any neighboring parities
66 : * (i.e. not the worst-case scenario) the array will be padded with zeros to
67 : * be the correct size. The first index may be 0 if the first component is
68 : * odd, but any following zeros are guaranteed to only have zeros following it.
69 : * The first returned `size_t` is the number of even components, while the
70 : * second is the number of odd components.
71 : */
72 : template <typename VariablesTags, bool IncludeZ = false>
73 : constexpr std::tuple<
74 : std::array<size_t,
75 : Variables<VariablesTags>::number_of_independent_components + 1>,
76 : size_t, size_t>
77 1 : compute_parity_list() {
78 : const size_t N =
79 : Variables<VariablesTags>::number_of_independent_components + 1;
80 : std::array<size_t, N> parity_run_lengths{};
81 :
82 : // x (spatial index 0 / spacetime index 1) always flips sign.
83 : // z (spatial index 2 / spacetime index 3) flips sign only when IncludeZ.
84 : const auto is_flipping_index = [](const IndexType index_type,
85 : const size_t index_value) {
86 : const bool is_x =
87 : (index_type == IndexType::Spacetime and index_value == 1) or
88 : (index_type == IndexType::Spatial and index_value == 0);
89 : const bool is_z =
90 : IncludeZ and
91 : ((index_type == IndexType::Spacetime and index_value == 3) or
92 : (index_type == IndexType::Spatial and index_value == 2));
93 : return is_x or is_z;
94 : };
95 :
96 : size_t run_index = 0;
97 : bool current_parity_is_even = true;
98 : tmpl::for_each<VariablesTags>([&parity_run_lengths, &run_index,
99 : ¤t_parity_is_even,
100 : &is_flipping_index]<typename TensorTag>(
101 : tmpl::type_<TensorTag> /*meta*/) {
102 : using tensor_type = typename TensorTag::type;
103 : constexpr auto index_types = tensor_type::index_types();
104 : constexpr size_t tensor_size = tensor_type::size();
105 :
106 : for (size_t component_index = 0; component_index < tensor_size;
107 : ++component_index) {
108 : const auto tensor_index = tensor_type::get_tensor_index(component_index);
109 : size_t flipping_index_count = 0;
110 : for (size_t index_position = 0; index_position < index_types.size();
111 : ++index_position) {
112 : if (is_flipping_index(gsl::at(index_types, index_position),
113 : gsl::at(tensor_index, index_position))) {
114 : ++flipping_index_count;
115 : }
116 : }
117 : // If current parity doesn't match last
118 : const bool component_is_even = (flipping_index_count % 2 == 0);
119 : if (component_is_even != current_parity_is_even) {
120 : ++run_index;
121 : current_parity_is_even = !current_parity_is_even;
122 : }
123 : gsl::at(parity_run_lengths, run_index) += 1;
124 : }
125 : });
126 :
127 : const auto [num_even, num_odd] = [&parity_run_lengths]<size_t... Is>(
128 : std::index_sequence<Is...> /*meta*/) {
129 : std::size_t even_count = ((Is % 2 == 0 ? parity_run_lengths[Is] : 0) + ...);
130 : std::size_t odd_count = ((Is % 2 == 1 ? parity_run_lengths[Is] : 0) + ...);
131 : return std::pair{even_count, odd_count};
132 : }(std::make_index_sequence<N>{});
133 :
134 : return {parity_run_lengths, num_even, num_odd};
135 : }
136 :
137 : /*!
138 : * \brief A compile-time function to determine parity of a tensor for an
139 : * axisymmetric spacetime.
140 : *
141 : * \see `compute_parity_list`
142 : */
143 : template <typename TensorType, bool IncludeZ = false>
144 : requires(tt::is_a_v<Tensor, TensorType>)
145 : constexpr std::tuple<std::array<size_t, TensorType::size() + 1>, size_t, size_t>
146 1 : compute_parity_list() {
147 : using tensor_type = Tensor<DataVector, typename TensorType::symmetry,
148 : typename TensorType::index_list>;
149 : using vars_list = ::Tags::convert_to_temp_tensors<tmpl::list<tensor_type>, 0>;
150 : return compute_parity_list<vars_list, IncludeZ>();
151 : }
152 : /*!
153 : * \brief Returns a compile-time array mapping each component index of
154 : * `TensorType` to its `Parity` (Even or Odd) for an axisymmetric spacetime.
155 : *
156 : * \see `compute_parity_list`
157 : */
158 : template <typename TensorType, bool IncludeZ = false>
159 : requires(tt::is_a_v<Tensor, TensorType>)
160 1 : constexpr std::array<Parity, TensorType::size()> make_component_parity_array() {
161 : constexpr auto parity_info = compute_parity_list<TensorType, IncludeZ>();
162 : constexpr auto parity_list = std::get<0>(parity_info);
163 : constexpr size_t N = TensorType::size();
164 : std::array<Parity, N> result{};
165 : size_t component = 0;
166 : bool is_even = true;
167 : for (size_t i = 0; component < N; ++i) {
168 : const size_t seg_size = parity_list[i];
169 : if (seg_size == 0) {
170 : if (is_even) {
171 : is_even = false;
172 : continue;
173 : } else {
174 : break;
175 : }
176 : }
177 : for (size_t k = 0; k < seg_size; ++k, ++component) {
178 : result[component] = is_even ? Parity::Even : Parity::Odd;
179 : }
180 : is_even = not is_even;
181 : }
182 : return result;
183 : }
184 : } // namespace Spectral
|
