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/DataVector.hpp"
10 : #include "DataStructures/Variables.hpp"
11 : #include "Evolution/DgSubcell/ReconstructionMethod.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 reconstruct_impl(gsl::span<double> dg_u,
31 : gsl::span<const double> subcell_u_times_projected_det_jac,
32 : const Mesh<Dim>& dg_mesh,
33 : const Index<Dim>& subcell_extents,
34 : ReconstructionMethod reconstruction_method,
35 : Spectral::Parity parity);
36 :
37 : /*!
38 : * \brief Reconstruct `subcell_u` onto the DG grid, sorting even- and
39 : * odd-parity components into separate batches when a ZernikeB1 basis is
40 : * present.
41 : *
42 : * For non-ZernikeB1 meshes this falls through to `reconstruct_impl` with
43 : * `Parity::Uninitialized`. The `TagList` must be the full Variables tag list
44 : * so that `Spectral::compute_parity_list` can determine per-component parity.
45 : */
46 : template <typename TagList, size_t Dim>
47 : void reconstruct_impl_with_tag_list(
48 : gsl::span<double> dg_u, gsl::span<const double> subcell_u,
49 : const Mesh<Dim>& dg_mesh, const Index<Dim>& subcell_extents,
50 : const ReconstructionMethod reconstruction_method) {
51 : if (dg_mesh.basis(0) == Spectral::Basis::ZernikeB1) {
52 : ASSERT(Variables<TagList>::number_of_independent_components *
53 : dg_mesh.number_of_grid_points() ==
54 : dg_u.size(),
55 : "Passed TagList does not have the same components, "
56 : << Variables<TagList>::number_of_independent_components
57 : << ", as dg_u holds, "
58 : << dg_u.size() / dg_mesh.number_of_grid_points());
59 : constexpr auto parity_info = Spectral::compute_parity_list<TagList>();
60 : constexpr auto parity_list = std::get<0>(parity_info);
61 : constexpr size_t num_even = std::get<1>(parity_info);
62 : constexpr size_t num_odd = std::get<2>(parity_info);
63 :
64 : const size_t num_dg_pts = dg_mesh.number_of_grid_points();
65 : const size_t num_subcell_pts = subcell_extents.product();
66 :
67 : // NOLINTNEXTLINE(modernize-avoid-c-arrays)
68 : auto buffer = cpp20::make_unique_for_overwrite<double[]>(
69 : (num_even + num_odd) * (num_subcell_pts + num_dg_pts));
70 : DataVector even_subcell_input{&buffer[0], num_even * num_subcell_pts};
71 : DataVector odd_subcell_input{&buffer[num_even * num_subcell_pts],
72 : num_odd * num_subcell_pts};
73 : DataVector even_dg_output{&buffer[(num_even + num_odd) * num_subcell_pts],
74 : num_even * num_dg_pts};
75 : DataVector odd_dg_output{
76 : &buffer[(num_even + num_odd) * num_subcell_pts + num_even * num_dg_pts],
77 : num_odd * num_dg_pts};
78 :
79 : // Sort subcell input into even/odd parity buffers
80 : const double* p_in = subcell_u.data();
81 : double* p_even_in = even_subcell_input.data();
82 : double* p_odd_in = odd_subcell_input.data();
83 : bool is_even = true;
84 : for (const size_t seg_size : parity_list) {
85 : if (seg_size == 0) {
86 : if (is_even) {
87 : is_even = false;
88 : continue;
89 : } else {
90 : break;
91 : }
92 : }
93 : if (is_even) {
94 : std::copy(p_in, p_in + seg_size * num_subcell_pts, // NOLINT
95 : p_even_in);
96 : p_even_in += seg_size * num_subcell_pts; // NOLINT
97 : } else {
98 : std::copy(p_in, p_in + seg_size * num_subcell_pts, // NOLINT
99 : p_odd_in);
100 : p_odd_in += seg_size * num_subcell_pts; // NOLINT
101 : }
102 : p_in += seg_size * num_subcell_pts; // NOLINT
103 : is_even = not is_even;
104 : }
105 :
106 : // Reconstruct each parity batch with the appropriate reconstruction
107 : // matrix
108 : if constexpr (num_even > 0) {
109 : reconstruct_impl(
110 : gsl::span<double>{even_dg_output.data(), even_dg_output.size()},
111 : gsl::span<const double>{even_subcell_input.data(),
112 : even_subcell_input.size()},
113 : dg_mesh, subcell_extents, reconstruction_method,
114 : Spectral::Parity::Even);
115 : }
116 : if constexpr (num_odd > 0) {
117 : reconstruct_impl(
118 : gsl::span<double>{odd_dg_output.data(), odd_dg_output.size()},
119 : gsl::span<const double>{odd_subcell_input.data(),
120 : odd_subcell_input.size()},
121 : dg_mesh, subcell_extents, reconstruction_method,
122 : Spectral::Parity::Odd);
123 : }
124 :
125 : // Reassemble output in original component order
126 : double* p_out = dg_u.data();
127 : const double* p_even_out = even_dg_output.data();
128 : const double* p_odd_out = odd_dg_output.data();
129 : is_even = true;
130 : for (const size_t seg_size : parity_list) {
131 : if (seg_size == 0) {
132 : if (is_even) {
133 : is_even = false;
134 : continue;
135 : } else {
136 : break;
137 : }
138 : }
139 : if (is_even) {
140 : // NOLINTNEXTLINE
141 : std::copy(p_even_out, p_even_out + seg_size * num_dg_pts, p_out);
142 : p_even_out += seg_size * num_dg_pts; // NOLINT
143 : } else {
144 : // NOLINTNEXTLINE
145 : std::copy(p_odd_out, p_odd_out + seg_size * num_dg_pts, p_out);
146 : p_odd_out += seg_size * num_dg_pts; // NOLINT
147 : }
148 : p_out += seg_size * num_dg_pts; // NOLINT
149 : is_even = not is_even;
150 : }
151 : return;
152 : }
153 : reconstruct_impl(dg_u, subcell_u, dg_mesh, subcell_extents,
154 : reconstruction_method, Spectral::Parity::Uninitialized);
155 : }
156 : } // namespace detail
157 :
158 : /// @{
159 : /*!
160 : * \ingroup DgSubcellGroup
161 : * \brief reconstruct the variable `subcell_u_times_projected_det_jac` onto the
162 : * DG grid `dg_mesh`.
163 : *
164 : * In general we wish that the reconstruction operator is the pseudo-inverse of
165 : * the projection operator. On curved meshes this means we either need to
166 : * compute a (time-dependent) reconstruction and projection matrix on each DG
167 : * element, or we expand the determinant of the Jacobian on the basis, accepting
168 : * the aliasing errors from that. We accept the aliasing errors in favor of the
169 : * significantly reduced computational overhead. This means that the projection
170 : * and reconstruction operators are only inverses of each other if both operate
171 : * on \f$u J\f$ where \f$u\f$ is the variable being projected and \f$J\f$ is the
172 : * determinant of the Jacobian. That is, the matrices are guaranteed to satisfy
173 : * \f$\mathcal{R}(\mathcal{P}(u J))=u J\f$. If the mesh is regular Cartesian,
174 : * then this isn't an issue. Furthermore, if we reconstruct
175 : * \f$uJ/\mathcal{P}(J)\f$ we again recover the exact DG solution. Doing the
176 : * latter has the advantage that, in general, we are ideally projecting to the
177 : * subcells much more often than reconstructing from them (a statement that we
178 : * would rather use DG more than the subcells).
179 : *
180 : * When the DG mesh uses a ZernikeB1 basis the Variables overloads deduce
181 : * per-component parity from the tag list automatically. The raw `DataVector`
182 : * overloads accepting a `Spectral::Parity` are for single-component data where
183 : * the caller already knows the parity. Only `DimByDim` reconstruction is
184 : * supported for ZernikeB1 meshes.
185 : */
186 : template <size_t Dim>
187 1 : DataVector reconstruct(
188 : const DataVector& subcell_u_times_projected_det_jac,
189 : const Mesh<Dim>& dg_mesh, const Index<Dim>& subcell_extents,
190 : ReconstructionMethod reconstruction_method,
191 : Spectral::Parity parity = Spectral::Parity::Uninitialized);
192 :
193 : template <size_t Dim>
194 1 : void reconstruct(gsl::not_null<DataVector*> dg_u,
195 : const DataVector& subcell_u_times_projected_det_jac,
196 : const Mesh<Dim>& dg_mesh, const Index<Dim>& subcell_extents,
197 : ReconstructionMethod reconstruction_method,
198 : Spectral::Parity parity = Spectral::Parity::Uninitialized);
199 :
200 : template <typename SubcellTagList, typename DgTagList, size_t Dim>
201 1 : void reconstruct(const gsl::not_null<Variables<DgTagList>*> dg_u,
202 : const Variables<SubcellTagList>& subcell_u,
203 : const Mesh<Dim>& dg_mesh, const Index<Dim>& subcell_extents,
204 : const ReconstructionMethod reconstruction_method) {
205 : ASSERT(subcell_u.number_of_grid_points() == subcell_extents.product(),
206 : "Incorrect subcell size of u: " << subcell_u.number_of_grid_points()
207 : << " but should be "
208 : << subcell_extents.product());
209 : if (UNLIKELY(dg_u->number_of_grid_points() !=
210 : dg_mesh.number_of_grid_points())) {
211 : dg_u->initialize(dg_mesh.number_of_grid_points(), 0.0);
212 : }
213 : detail::reconstruct_impl_with_tag_list<DgTagList>(
214 : gsl::span<double>{dg_u->data(), dg_u->size()},
215 : gsl::span<const double>{subcell_u.data(), subcell_u.size()}, dg_mesh,
216 : subcell_extents, reconstruction_method);
217 : }
218 :
219 : template <typename TagList, size_t Dim>
220 1 : Variables<TagList> reconstruct(
221 : const Variables<TagList>& subcell_u, const Mesh<Dim>& dg_mesh,
222 : const Index<Dim>& subcell_extents,
223 : const ReconstructionMethod reconstruction_method) {
224 : Variables<TagList> dg_u(dg_mesh.number_of_grid_points());
225 : reconstruct(make_not_null(&dg_u), subcell_u, dg_mesh, subcell_extents,
226 : reconstruction_method);
227 : return dg_u;
228 : }
229 : /// @}
230 : } // namespace evolution::dg::subcell::fd
|
