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