Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <algorithm>
7 : #include <cstddef>
8 :
9 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
10 : #include "DataStructures/DataVector.hpp"
11 : #include "DataStructures/Variables.hpp"
12 : #include "NumericalAlgorithms/Spectral/Basis.hpp"
13 : #include "NumericalAlgorithms/Spectral/Parity.hpp"
14 : #include "NumericalAlgorithms/Spectral/ParityFromSymmetry.hpp"
15 : #include "Utilities/ErrorHandling/Assert.hpp"
16 : #include "Utilities/Gsl.hpp"
17 : #include "Utilities/MemoryHelpers.hpp"
18 : #include "Utilities/TMPL.hpp"
19 :
20 : /// \cond
21 : template <size_t>
22 : class Index;
23 : template <size_t>
24 : class Mesh;
25 : /// \endcond
26 :
27 : namespace evolution::dg::subcell::fd {
28 : namespace detail {
29 : template <size_t Dim>
30 : void project_impl(gsl::span<double> subcell_u, gsl::span<const double> dg_u,
31 : const Mesh<Dim>& dg_mesh, const Index<Dim>& subcell_extents,
32 : Spectral::Parity parity);
33 : template <size_t Dim>
34 : void project_to_faces_impl(gsl::span<double> subcell_u,
35 : gsl::span<const double> dg_u,
36 : const Mesh<Dim>& dg_mesh,
37 : const Index<Dim>& subcell_extents,
38 : const size_t& face_direction,
39 : Spectral::Parity parity);
40 :
41 : /*!
42 : * \brief Project `dg_u` onto the subcell grid, sorting even- and odd-parity
43 : * components into separate batches when a ZernikeB1 basis is present.
44 : *
45 : * For non-ZernikeB1 meshes this falls through to the default implementation.
46 : * The `TagList` must be the full Variables tag list.
47 : */
48 : template <typename TagList, size_t Dim>
49 : void project_impl_with_tag_list(gsl::span<double> subcell_u,
50 : gsl::span<const double> dg_u,
51 : const Mesh<Dim>& dg_mesh,
52 : const Index<Dim>& subcell_extents) {
53 : if (dg_mesh.basis(0) == Spectral::Basis::ZernikeB1) {
54 : ASSERT(Variables<TagList>::number_of_independent_components *
55 : dg_mesh.number_of_grid_points() ==
56 : dg_u.size(),
57 : "Passed TagList does not have the same components, "
58 : << Variables<TagList>::number_of_independent_components
59 : << ", as dg_u holds, "
60 : << dg_u.size() / dg_mesh.number_of_grid_points());
61 : constexpr auto parity_info = Spectral::compute_parity_list<TagList>();
62 : constexpr auto parity_list = std::get<0>(parity_info);
63 : constexpr size_t num_even = std::get<1>(parity_info);
64 : constexpr size_t num_odd = std::get<2>(parity_info);
65 :
66 : const size_t num_dg_pts = dg_mesh.number_of_grid_points();
67 : const size_t num_subcell_pts = subcell_extents.product();
68 :
69 : // NOLINTNEXTLINE(modernize-avoid-c-arrays)
70 : auto buffer = cpp20::make_unique_for_overwrite<double[]>(
71 : (num_even + num_odd) * (num_dg_pts + num_subcell_pts));
72 : DataVector even_input{&buffer[0], num_even * num_dg_pts};
73 : DataVector odd_input(&buffer[num_even * num_dg_pts], num_odd * num_dg_pts);
74 : DataVector even_output(&buffer[(num_even + num_odd) * num_dg_pts],
75 : num_even * num_subcell_pts);
76 : DataVector odd_output(
77 : &buffer[(num_even + num_odd) * num_dg_pts + num_even * num_subcell_pts],
78 : num_odd * num_subcell_pts);
79 :
80 : // Sort input components into even/odd parity buffers
81 : const double* p_in = dg_u.data();
82 : double* p_even_in = even_input.data();
83 : double* p_odd_in = odd_input.data();
84 : bool is_even = true;
85 : for (const size_t seg_size : parity_list) {
86 : if (seg_size == 0) {
87 : if (is_even) {
88 : is_even = false;
89 : continue;
90 : } else {
91 : break;
92 : }
93 : }
94 : if (is_even) {
95 : std::copy(p_in, p_in + seg_size * num_dg_pts, p_even_in); // NOLINT
96 : p_even_in += seg_size * num_dg_pts; // NOLINT
97 : } else {
98 : std::copy(p_in, p_in + seg_size * num_dg_pts, p_odd_in); // NOLINT
99 : p_odd_in += seg_size * num_dg_pts; // NOLINT
100 : }
101 : p_in += seg_size * num_dg_pts; // NOLINT
102 : is_even = not is_even;
103 : }
104 :
105 : // Project each parity batch with the appropriate projection matrix
106 : if constexpr (num_even > 0) {
107 : project_impl(
108 : gsl::span<double>{even_output.data(), even_output.size()},
109 : gsl::span<const double>{even_input.data(), even_input.size()},
110 : dg_mesh, subcell_extents, Spectral::Parity::Even);
111 : }
112 : if constexpr (num_odd > 0) {
113 : project_impl(gsl::span<double>{odd_output.data(), odd_output.size()},
114 : gsl::span<const double>{odd_input.data(), odd_input.size()},
115 : dg_mesh, subcell_extents, Spectral::Parity::Odd);
116 : }
117 :
118 : // Reassemble output in original component order
119 : double* p_out = subcell_u.data();
120 : const double* p_even_out = even_output.data();
121 : const double* p_odd_out = odd_output.data();
122 : is_even = true;
123 : for (const size_t seg_size : parity_list) {
124 : if (seg_size == 0) {
125 : if (is_even) {
126 : is_even = false;
127 : continue;
128 : } else {
129 : break;
130 : }
131 : }
132 : if (is_even) {
133 : // NOLINTNEXTLINE
134 : std::copy(p_even_out, p_even_out + seg_size * num_subcell_pts, p_out);
135 : p_even_out += seg_size * num_subcell_pts; // NOLINT
136 : } else {
137 : // NOLINTNEXTLINE
138 : std::copy(p_odd_out, p_odd_out + seg_size * num_subcell_pts, p_out);
139 : p_odd_out += seg_size * num_subcell_pts; // NOLINT
140 : }
141 : p_out += seg_size * num_subcell_pts; // NOLINT
142 : is_even = not is_even;
143 : }
144 : return;
145 : }
146 : project_impl(subcell_u, dg_u, dg_mesh, subcell_extents,
147 : Spectral::Parity::Uninitialized);
148 : }
149 : } // namespace detail
150 :
151 : /// @{
152 : /*!
153 : * \ingroup DgSubcellGroup
154 : * \brief Project the variable `dg_u` onto the subcell grid with extents
155 : * `subcell_extents`.
156 : *
157 : * When the DG mesh uses a ZernikeB1 basis the Variables overloads deduce
158 : * per-component parity from the tag list automatically. The raw `DataVector`
159 : * overloads accepting a `Spectral::Parity` are for single-component data
160 : * where the caller already knows the parity; the overloads accepting a
161 : * `tmpl::list<TagList>` meta-parameter project a multi-component DataVector
162 : * whose tensor-parity structure is encoded in `TagList`.
163 : *
164 : * \note In the return-by-`gsl::not_null` with `Variables` interface, the
165 : * `SubcellTagList` and the `DgTagList` must be the same when all tag prefixes
166 : * are removed.
167 : */
168 : template <size_t Dim>
169 1 : DataVector project(const DataVector& dg_u, const Mesh<Dim>& dg_mesh,
170 : const Index<Dim>& subcell_extents,
171 : Spectral::Parity parity = Spectral::Parity::Uninitialized);
172 :
173 : template <size_t Dim>
174 1 : void project(gsl::not_null<DataVector*> subcell_u, const DataVector& dg_u,
175 : const Mesh<Dim>& dg_mesh, const Index<Dim>& subcell_extents,
176 : Spectral::Parity parity = Spectral::Parity::Uninitialized);
177 :
178 : template <typename TagList, size_t Dim>
179 1 : void project(const gsl::not_null<DataVector*> subcell_u, const DataVector& dg_u,
180 : const Mesh<Dim>& dg_mesh, const Index<Dim>& subcell_extents,
181 : TagList /*meta*/) {
182 : ASSERT(dg_u.size() % dg_mesh.number_of_grid_points() == 0,
183 : "The vector dg_u must have size that is a multiple of the number of "
184 : "grid points "
185 : << dg_mesh.number_of_grid_points() << " but got " << dg_u.size());
186 : subcell_u->destructive_resize(subcell_extents.product() * dg_u.size() /
187 : dg_mesh.number_of_grid_points());
188 : detail::project_impl_with_tag_list<TagList>(
189 : gsl::span<double>{subcell_u->data(), subcell_u->size()},
190 : gsl::span<const double>{dg_u.data(), dg_u.size()}, dg_mesh,
191 : subcell_extents);
192 : }
193 :
194 : template <typename TagList, size_t Dim>
195 1 : DataVector project(const DataVector& dg_u, const Mesh<Dim>& dg_mesh,
196 : const Index<Dim>& subcell_extents, TagList /*meta*/) {
197 : ASSERT(dg_u.size() % dg_mesh.number_of_grid_points() == 0,
198 : "The vector dg_u must have size that is a multiple of the number of "
199 : "grid points "
200 : << dg_mesh.number_of_grid_points() << " but got " << dg_u.size());
201 : DataVector subcell_u(subcell_extents.product() * dg_u.size() /
202 : dg_mesh.number_of_grid_points());
203 : project(make_not_null(&subcell_u), dg_u, dg_mesh, subcell_extents, TagList{});
204 : return subcell_u;
205 : }
206 :
207 : template <typename SubcellTagList, typename DgTagList, size_t Dim>
208 1 : void project(const gsl::not_null<Variables<SubcellTagList>*> subcell_u,
209 : const Variables<DgTagList>& dg_u, const Mesh<Dim>& dg_mesh,
210 : const Index<Dim>& subcell_extents) {
211 : static_assert(
212 : std::is_same_v<
213 : tmpl::transform<SubcellTagList,
214 : tmpl::bind<db::remove_all_prefixes, tmpl::_1>>,
215 : tmpl::transform<DgTagList,
216 : tmpl::bind<db::remove_all_prefixes, tmpl::_1>>>,
217 : "DG and subcell tag lists must be the same once prefix tags "
218 : "are removed.");
219 : ASSERT(dg_u.number_of_grid_points() == dg_mesh.number_of_grid_points(),
220 : "dg_u has incorrect size " << dg_u.number_of_grid_points()
221 : << " since the mesh is size "
222 : << dg_mesh.number_of_grid_points());
223 : if (UNLIKELY(subcell_u->number_of_grid_points() !=
224 : subcell_extents.product())) {
225 : subcell_u->initialize(subcell_extents.product());
226 : }
227 : detail::project_impl_with_tag_list<DgTagList>(
228 : gsl::span<double>{subcell_u->data(), subcell_u->size()},
229 : gsl::span<const double>{dg_u.data(), dg_u.size()}, dg_mesh,
230 : subcell_extents);
231 : }
232 :
233 : template <typename TagList, size_t Dim>
234 1 : Variables<TagList> project(const Variables<TagList>& dg_u,
235 : const Mesh<Dim>& dg_mesh,
236 : const Index<Dim>& subcell_extents) {
237 : Variables<TagList> subcell_u(subcell_extents.product());
238 : project(make_not_null(&subcell_u), dg_u, dg_mesh, subcell_extents);
239 : return subcell_u;
240 : }
241 :
242 : template <size_t Dim>
243 1 : DataVector project_to_faces(const DataVector& dg_u, const Mesh<Dim>& dg_mesh,
244 : const Index<Dim>& subcell_extents,
245 : const size_t& face_direction,
246 : Spectral::Parity parity);
247 :
248 : template <size_t Dim>
249 1 : void project_to_faces(gsl::not_null<DataVector*> subcell_u,
250 : const DataVector& dg_u, const Mesh<Dim>& dg_mesh,
251 : const Index<Dim>& subcell_extents,
252 : const size_t& face_direction, Spectral::Parity parity);
253 :
254 : template <typename SubcellTagList, typename DgTagList, size_t Dim>
255 1 : void project_to_faces(const gsl::not_null<Variables<SubcellTagList>*> subcell_u,
256 : const Variables<DgTagList>& dg_u,
257 : const Mesh<Dim>& dg_mesh,
258 : const Index<Dim>& subcell_extents,
259 : const size_t& face_direction, Spectral::Parity parity) {
260 : static_assert(
261 : std::is_same_v<
262 : tmpl::transform<SubcellTagList,
263 : tmpl::bind<db::remove_all_prefixes, tmpl::_1>>,
264 : tmpl::transform<DgTagList,
265 : tmpl::bind<db::remove_all_prefixes, tmpl::_1>>>,
266 : "DG and subcell tag lists must be the same once prefix tags "
267 : "are removed.");
268 : ASSERT(dg_u.number_of_grid_points() == dg_mesh.number_of_grid_points(),
269 : "dg_u has incorrect size " << dg_u.number_of_grid_points()
270 : << " since the mesh is size "
271 : << dg_mesh.number_of_grid_points());
272 : if (UNLIKELY(subcell_u->number_of_grid_points() !=
273 : subcell_extents.product())) {
274 : subcell_u->initialize(subcell_extents.product());
275 : }
276 : detail::project_to_faces_impl(
277 : gsl::span<double>{subcell_u->data(), subcell_u->size()},
278 : gsl::span<const double>{dg_u.data(), dg_u.size()}, dg_mesh,
279 : subcell_extents, face_direction, parity);
280 : }
281 :
282 : template <typename TagList, size_t Dim>
283 1 : Variables<TagList> project_to_faces(const Variables<TagList>& dg_u,
284 : const Mesh<Dim>& dg_mesh,
285 : const Index<Dim>& subcell_extents,
286 : const size_t& face_direction,
287 : Spectral::Parity parity) {
288 : Variables<TagList> subcell_u(subcell_extents.product());
289 : project_to_faces(make_not_null(&subcell_u), dg_u, dg_mesh, subcell_extents,
290 : face_direction, parity);
291 : return subcell_u;
292 : }
293 : /// @}
294 : } // namespace evolution::dg::subcell::fd
|
