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 <optional>
9 : #include <unordered_set>
10 : #include <utility>
11 : #include <vector>
12 :
13 : #include "DataStructures/DataBox/DataBox.hpp"
14 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
15 : #include "DataStructures/DataBox/Prefixes.hpp"
16 : #include "DataStructures/DataVector.hpp"
17 : #include "DataStructures/Index.hpp"
18 : #include "Domain/Structure/Direction.hpp"
19 : #include "Domain/Structure/DirectionMap.hpp"
20 : #include "Domain/Structure/Element.hpp"
21 : #include "Domain/Structure/ElementId.hpp"
22 : #include "Domain/Structure/OrientationMapHelpers.hpp"
23 : #include "Domain/Tags.hpp"
24 : #include "Domain/Tags/NeighborMesh.hpp"
25 : #include "Evolution/DgSubcell/Mesh.hpp"
26 : #include "Evolution/DgSubcell/Projection.hpp"
27 : #include "Evolution/DgSubcell/RdmpTci.hpp"
28 : #include "Evolution/DgSubcell/RdmpTciData.hpp"
29 : #include "Evolution/DgSubcell/SliceData.hpp"
30 : #include "Evolution/DgSubcell/SubcellOptions.hpp"
31 : #include "Evolution/DgSubcell/Tags/DataForRdmpTci.hpp"
32 : #include "Evolution/DgSubcell/Tags/Interpolators.hpp"
33 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
34 : #include "Evolution/DgSubcell/Tags/SubcellOptions.hpp"
35 : #include "Evolution/DiscontinuousGalerkin/OnlyDgBlockIds.hpp"
36 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
37 : #include "Utilities/Algorithm.hpp"
38 : #include "Utilities/ErrorHandling/Assert.hpp"
39 : #include "Utilities/Gsl.hpp"
40 : #include "Utilities/TMPL.hpp"
41 :
42 : namespace evolution::dg::subcell {
43 : /*!
44 : * \brief Add local data for our and our neighbor's relaxed discrete maximum
45 : * principle troubled-cell indicator, and for reconstruction.
46 : *
47 : * The local maximum and minimum of the evolved variables is added to
48 : * `Tags::DataForRdmpTci` for the RDMP TCI. Then the
49 : * data needed by neighbor elements to do reconstruction on the FD grid is sent.
50 : * The data to be sent is computed in the mutator
51 : * `Metavariables::SubcellOptions::GhostVariables`, which returns a
52 : * `Variables` of the tensors to send to the neighbors. The main reason for
53 : * having the mutator `GhostVariables` is to allow sending primitive or
54 : * characteristic variables for reconstruction.
55 : *
56 : * \note If all neighbors are using DG (i.e. none of their meshes is an FD mesh)
57 : * then we send our DG volume data _without_ orienting it. This elides the
58 : * expense of projection and slicing. If any
59 : * neighbors are doing FD, we project and slice to all neighbors. A future
60 : * optimization would be to measure the cost of slicing data, and figure out how
61 : * many neighbors need to be doing FD before it's worth projecting and slicing
62 : * vs. just projecting to the ghost cells. Another optimization is to always
63 : * send DG volume data to neighbors that we think are doing DG, so that we can
64 : * elide any slicing or projection cost in that direction.
65 : *
66 : * No data is prepared for neighbors that can never use subcell. If our own
67 : * block is DG-only (see `evolution::dg::Tags::OnlyDgBlockIds`), all of our
68 : * neighbors border a DG-only block and so never use subcell. In this case
69 : * data is never prepared and `all_neighbor_data_for_reconstruction` is left
70 : * empty. Otherwise, we skip the directions in which all neighbors are in
71 : * DG-only blocks, meaing we prepare only for directions in which there are
72 : * neighbors that could use subcell (this check currently still sends extra
73 : * data to DG-only elements buffering subcell and DG-only blocks as
74 : * `OnlyDgBlockIds` does not have element granulaity).
75 : */
76 : template <typename Metavariables, typename DbTagsList, size_t Dim>
77 1 : void prepare_neighbor_data(
78 : const gsl::not_null<DirectionMap<Dim, DataVector>*>
79 : all_neighbor_data_for_reconstruction,
80 : const gsl::not_null<std::optional<Mesh<Dim>>*> ghost_data_mesh,
81 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
82 : [[maybe_unused]] const Variables<db::wrap_tags_in<
83 : ::Tags::Flux, typename Metavariables::system::flux_variables,
84 : tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes) {
85 : const Mesh<Dim>& dg_mesh = db::get<::domain::Tags::Mesh<Dim>>(*box);
86 : const Mesh<Dim>& subcell_mesh =
87 : db::get<::evolution::dg::subcell::Tags::Mesh<Dim>>(*box);
88 : const auto& element = db::get<::domain::Tags::Element<Dim>>(*box);
89 : const std::vector<size_t>& only_dg_block_ids =
90 : db::get<evolution::dg::Tags::OnlyDgBlockIds<Dim>>(*box);
91 :
92 : all_neighbor_data_for_reconstruction->clear();
93 : *ghost_data_mesh = dg_mesh;
94 : // Ghost data is only used to fill the ghost cells of neighbors doing FD, and
95 : // elements in or bordering a DG-only block never use subcell (see
96 : // `subcell::Actions::Initialize`). So if our block is DG-only none of our
97 : // neighbors needs ghost data, and otherwise we skip the directions in which
98 : // all neighbors are in DG-only blocks.
99 : if (alg::found(only_dg_block_ids, element.id().block_id())) {
100 : return;
101 : }
102 : std::unordered_set<Direction<Dim>> directions_to_slice{};
103 : for (const auto& [direction, neighbors] : element.neighbors()) {
104 : if (not alg::all_of(neighbors.ids(), [&only_dg_block_ids](
105 : const ElementId<Dim>& neighbor) {
106 : return alg::found(only_dg_block_ids, neighbor.block_id());
107 : })) {
108 : directions_to_slice.insert(direction);
109 : }
110 : }
111 : if (directions_to_slice.empty()) {
112 : return;
113 : }
114 :
115 : const auto& neighbor_meshes = db::get<domain::Tags::NeighborMesh<Dim>>(*box);
116 : const subcell::RdmpTciData& rdmp_data =
117 : db::get<subcell::Tags::DataForRdmpTci>(*box);
118 : const ::fd::DerivativeOrder fd_derivative_order =
119 : get<evolution::dg::subcell::Tags::SubcellOptions<Dim>>(*box)
120 : .finite_difference_derivative_order();
121 : const size_t rdmp_size = rdmp_data.max_variables_values.size() +
122 : rdmp_data.min_variables_values.size();
123 : const size_t extra_size_for_ghost_data =
124 : (fd_derivative_order == ::fd::DerivativeOrder::Two
125 : ? 0
126 : : volume_fluxes.size());
127 :
128 : if (DataVector ghost_variables =
129 : [&box, &extra_size_for_ghost_data, &rdmp_size, &volume_fluxes]() {
130 : if (extra_size_for_ghost_data == 0) {
131 : return db::mutate_apply(
132 : typename Metavariables::SubcellOptions::GhostVariables{}, box,
133 : rdmp_size);
134 : } else {
135 : DataVector ghost_vars = db::mutate_apply(
136 : typename Metavariables::SubcellOptions::GhostVariables{}, box,
137 : rdmp_size + extra_size_for_ghost_data);
138 : std::copy(
139 : volume_fluxes.data(),
140 : std::next(volume_fluxes.data(),
141 : static_cast<std::ptrdiff_t>(volume_fluxes.size())),
142 : std::prev(ghost_vars.end(),
143 : static_cast<std::ptrdiff_t>(
144 : rdmp_size + extra_size_for_ghost_data)));
145 : return ghost_vars;
146 : }
147 : }();
148 : (not fd::dg_mesh_supports_subcell(dg_mesh) or
149 : alg::all_of(neighbor_meshes,
150 : [](const auto& directional_element_id_and_mesh) {
151 : ASSERT(directional_element_id_and_mesh.second.basis(0) !=
152 : Spectral::Basis::Chebyshev,
153 : "Don't yet support Chebyshev basis with DG-FD");
154 : return directional_element_id_and_mesh.second.basis(0) !=
155 : Spectral::Basis::FiniteDifference;
156 : }))) {
157 : *ghost_data_mesh = dg_mesh;
158 : const size_t total_to_slice = directions_to_slice.size();
159 : size_t slice_count = 0;
160 : for (const auto& direction : directions_to_slice) {
161 : ++slice_count;
162 : // Move instead of copy on the last iteration. Elides a memory
163 : // allocation and copy.
164 : [[maybe_unused]] const auto insert_result =
165 : UNLIKELY(slice_count == total_to_slice)
166 : ? all_neighbor_data_for_reconstruction->insert_or_assign(
167 : direction, std::move(ghost_variables))
168 : : all_neighbor_data_for_reconstruction->insert_or_assign(
169 : direction, ghost_variables);
170 : ASSERT(all_neighbor_data_for_reconstruction->size() == total_to_slice or
171 : insert_result.second,
172 : "Failed to insert the neighbor data in direction " << direction);
173 : }
174 : } else {
175 : *ghost_data_mesh = subcell_mesh;
176 : const size_t ghost_zone_size =
177 : Metavariables::SubcellOptions::ghost_zone_size(*box);
178 :
179 : const DataVector data_to_project{};
180 : make_const_view(make_not_null(&data_to_project), ghost_variables, 0,
181 : ghost_variables.size() - rdmp_size);
182 : const DataVector projected_data = evolution::dg::subcell::fd::project(
183 : data_to_project, dg_mesh, subcell_mesh.extents(),
184 : typename Metavariables::SubcellOptions::GhostVariables::
185 : ghost_variables_tag_list{});
186 : typename evolution::dg::subcell::Tags::InterpolatorsFromFdToNeighborFd<
187 : Dim>::type directions_not_to_slice{};
188 : for (const auto& [directional_element_id, _] :
189 : db::get<evolution::dg::subcell::Tags::InterpolatorsFromFdToNeighborFd<
190 : Dim>>(*box)) {
191 : directions_not_to_slice[directional_element_id] = std::nullopt;
192 : }
193 : *all_neighbor_data_for_reconstruction = subcell::slice_data(
194 : projected_data, db::get<subcell::Tags::Mesh<Dim>>(*box).extents(),
195 : ghost_zone_size, directions_to_slice, rdmp_size,
196 : directions_not_to_slice);
197 : for (const auto& [directional_element_id, interpolator] :
198 : db::get<evolution::dg::subcell::Tags::InterpolatorsFromDgToNeighborFd<
199 : Dim>>(*box)) {
200 : // NOTE: no orienting is needed, that's done by the interpolation.
201 : ASSERT(interpolator.has_value(),
202 : "All interpolators must have values. The std::optional is used "
203 : "for some local optimizations.");
204 : if (not directions_to_slice.contains(
205 : directional_element_id.direction())) {
206 : continue;
207 : }
208 : DataVector& data_in_direction = all_neighbor_data_for_reconstruction->at(
209 : directional_element_id.direction());
210 : gsl::span result_in_direction{data_in_direction.data(),
211 : data_in_direction.size() - rdmp_size};
212 : interpolator->interpolate(
213 : make_not_null(&result_in_direction),
214 : {data_to_project.data(), data_to_project.size()});
215 : }
216 : }
217 :
218 : // Note: The RDMP TCI data must be filled by the TCIs before getting to this
219 : // call. That means once in the initial data and then in both the DG and FD
220 : // TCIs.
221 : for (auto& [direction, neighbor_data] :
222 : *all_neighbor_data_for_reconstruction) {
223 : // Add the RDMP TCI data to what we will be sending.
224 : // Note that this is added _after_ the reconstruction data has been
225 : // re-oriented (in the case where we have neighbors doing FD).
226 : std::copy(rdmp_data.max_variables_values.begin(),
227 : rdmp_data.max_variables_values.end(),
228 : std::prev(neighbor_data.end(), static_cast<int>(rdmp_size)));
229 : std::copy(
230 : rdmp_data.min_variables_values.begin(),
231 : rdmp_data.min_variables_values.end(),
232 : std::prev(neighbor_data.end(),
233 : static_cast<int>(rdmp_data.min_variables_values.size())));
234 : }
235 : }
236 : } // namespace evolution::dg::subcell
|
