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 : #include <deque>
9 : #include <iterator>
10 : #include <optional>
11 : #include <tuple>
12 : #include <type_traits>
13 : #include <utility>
14 :
15 : #include "DataStructures/DataBox/DataBox.hpp"
16 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
17 : #include "DataStructures/DataBox/Prefixes.hpp"
18 : #include "DataStructures/DataVector.hpp"
19 : #include "DataStructures/TaggedTuple.hpp"
20 : #include "Domain/Structure/Direction.hpp"
21 : #include "Domain/Structure/DirectionalId.hpp"
22 : #include "Domain/Structure/DirectionalIdMap.hpp"
23 : #include "Domain/Structure/Element.hpp"
24 : #include "Domain/Structure/ElementId.hpp"
25 : #include "Domain/Tags.hpp"
26 : #include "Evolution/DgSubcell/Actions/Labels.hpp"
27 : #include "Evolution/DgSubcell/ActiveGrid.hpp"
28 : #include "Evolution/DgSubcell/GhostData.hpp"
29 : #include "Evolution/DgSubcell/Mesh.hpp"
30 : #include "Evolution/DgSubcell/NeighborRdmpAndVolumeData.hpp"
31 : #include "Evolution/DgSubcell/Projection.hpp"
32 : #include "Evolution/DgSubcell/RdmpTci.hpp"
33 : #include "Evolution/DgSubcell/RdmpTciData.hpp"
34 : #include "Evolution/DgSubcell/SubcellOptions.hpp"
35 : #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
36 : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
37 : #include "Evolution/DgSubcell/Tags/DataForRdmpTci.hpp"
38 : #include "Evolution/DgSubcell/Tags/DidRollback.hpp"
39 : #include "Evolution/DgSubcell/Tags/GhostDataForReconstruction.hpp"
40 : #include "Evolution/DgSubcell/Tags/Interpolators.hpp"
41 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
42 : #include "Evolution/DgSubcell/Tags/MeshForGhostData.hpp"
43 : #include "Evolution/DgSubcell/Tags/SubcellOptions.hpp"
44 : #include "Evolution/DgSubcell/Tags/TciStatus.hpp"
45 : #include "Evolution/DiscontinuousGalerkin/InboxTags.hpp"
46 : #include "Evolution/DiscontinuousGalerkin/OnlyDgBlockIds.hpp"
47 : #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp"
48 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
49 : #include "Parallel/AlgorithmExecution.hpp"
50 : #include "Parallel/GlobalCache.hpp"
51 : #include "ParallelAlgorithms/Actions/Goto.hpp"
52 : #include "Time/Actions/SelfStartActions.hpp"
53 : #include "Time/History.hpp"
54 : #include "Utilities/Algorithm.hpp"
55 : #include "Utilities/ContainerHelpers.hpp"
56 : #include "Utilities/ErrorHandling/Assert.hpp"
57 : #include "Utilities/TMPL.hpp"
58 :
59 : /// \cond
60 : namespace Tags {
61 : template <typename Tag>
62 : struct HistoryEvolvedVariables;
63 : struct TimeStepId;
64 : } // namespace Tags
65 : /// \endcond
66 :
67 : namespace evolution::dg::subcell::Actions {
68 : /*!
69 : * \brief Run the troubled-cell indicator on the candidate solution and perform
70 : * the time step rollback if needed.
71 : *
72 : * Elements that cannot use subcell (e.g. non-hypercube topology or DG-only
73 : * blocks) skip the TCI entirely, set the TCI decision to 0, clear ghost data,
74 : * and continue.
75 : *
76 : * Interior cells are marked as troubled if
77 : * `subcell_options.always_use_subcells()` is `true`, or if either the RDMP
78 : * troubled-cell indicator (TCI) or the TciMutator reports the cell is
79 : * troubled. Exterior cells are marked as troubled only if
80 : * `Metavariables::SubcellOptions::subcell_enabled_at_external_boundary` is
81 : * `true`.
82 : *
83 : * The troubled-cell indicator (TCI) given by the mutator `TciMutator` can
84 : * mutate tags in the DataBox, but should do so cautiously. The main reason that
85 : * this is a mutator is because primitive variables, such as the pressure, are
86 : * used to check if the solution is physical. In the relativistic case, even
87 : * just whether or not the primitive variables can be recovered is used as a
88 : * condition. Note that the evolved variables are projected to the subcells
89 : * _after_ the TCI is called and marks the cell as troubled.
90 : *
91 : * After rollback, the subcell scheme must project the DG boundary corrections
92 : * \f$G\f$ to the subcells for the scheme to be conservative. The subcell
93 : * actions know if a rollback was done because the local mortar data would
94 : * already be computed.
95 : */
96 : template <typename TciMutator>
97 1 : struct TciAndRollback {
98 : template <typename DbTags, typename... InboxTags, typename Metavariables,
99 : typename ArrayIndex, typename ActionList,
100 : typename ParallelComponent, size_t Dim = Metavariables::volume_dim>
101 0 : static Parallel::iterable_action_return_t apply(
102 : db::DataBox<DbTags>& box,
103 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
104 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
105 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
106 : const ParallelComponent* const /*meta*/) {
107 : static_assert(
108 : tmpl::count_if<
109 : ActionList,
110 : std::is_same<tmpl::_1, tmpl::pin<TciAndRollback>>>::value == 1,
111 : "Must have the TciAndRollback action exactly once in the action list "
112 : "of a phase.");
113 : static_assert(
114 : tmpl::count_if<
115 : ActionList,
116 : std::is_same<tmpl::_1,
117 : tmpl::pin<::Actions::Label<
118 : evolution::dg::subcell::Actions::Labels::
119 : BeginSubcellAfterDgRollback>>>>::value == 1,
120 : "Must have the BeginSubcellAfterDgRollback label exactly once in the "
121 : "action list of a phase.");
122 :
123 : using variables_tag = typename Metavariables::system::variables_tag;
124 :
125 : const ActiveGrid active_grid = db::get<Tags::ActiveGrid>(box);
126 : ASSERT(active_grid == ActiveGrid::Dg,
127 : "Must be using DG when calling TciAndRollback action.");
128 :
129 : const Element<Dim>& element = db::get<::domain::Tags::Element<Dim>>(box);
130 : const bool cell_has_external_boundary =
131 : not element.external_boundaries().empty();
132 :
133 : constexpr bool subcell_enabled_at_external_boundary =
134 : Metavariables::SubcellOptions::subcell_enabled_at_external_boundary;
135 : const Mesh<Dim>& dg_mesh = db::get<::domain::Tags::Mesh<Dim>>(box);
136 : const Mesh<Dim>& subcell_mesh = db::get<Tags::Mesh<Dim>>(box);
137 :
138 : const SubcellOptions& subcell_options =
139 : db::get<Tags::SubcellOptions<Dim>>(box);
140 : const std::vector<size_t>& only_dg_block_ids =
141 : db::get<evolution::dg::Tags::OnlyDgBlockIds<Dim>>(box);
142 : bool cell_is_troubled =
143 : subcell_options.always_use_subcells() or
144 : (subcell_options.use_halo() and [&box]() -> bool {
145 : for (const auto& [_, neighbor_decision] :
146 : db::get<evolution::dg::subcell::Tags::NeighborTciDecisions<Dim>>(
147 : box)) {
148 : if (neighbor_decision != 0) {
149 : return true;
150 : }
151 : }
152 : return false;
153 : }());
154 :
155 : // Loop over block neighbors and if neighbor id is inside of
156 : // only_dg_block_ids, then bordering DG-only block
157 : const bool bordering_dg_block = alg::any_of(
158 : element.neighbors(),
159 : [&only_dg_block_ids](const auto& direction_and_neighbor) {
160 : const size_t first_block_id =
161 : direction_and_neighbor.second.ids().begin()->block_id();
162 : return alg::found(only_dg_block_ids, first_block_id);
163 : });
164 :
165 : // Subcell is allowed in the element if 3 conditions are met:
166 : // (i) The DG mesh topology supports subcell
167 : // (ii) The current element block id is not marked as DG only
168 : // (iii) The current element is not bordering a DG only block.
169 : const bool subcell_allowed_in_element =
170 : fd::dg_mesh_supports_subcell(dg_mesh) and
171 : not alg::found(only_dg_block_ids, element.id().block_id()) and
172 : not bordering_dg_block;
173 :
174 : // Elements that can never use subcell (e.g. non-hypercube topology or
175 : // DG-only blocks) should skip the TCI entirely, since the TCI projects
176 : // to the subcell mesh which is unsupported for these elements.
177 : if (not subcell_allowed_in_element) {
178 : db::mutate<Tags::TciDecision,
179 : subcell::Tags::GhostDataForReconstruction<Dim>>(
180 : [](const gsl::not_null<int*> tci_decision_ptr,
181 : const auto neighbor_data_ptr) {
182 : *tci_decision_ptr = 0;
183 : neighbor_data_ptr->clear();
184 : },
185 : make_not_null(&box));
186 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
187 : }
188 :
189 : // The reason we pass in the persson_exponent explicitly instead of
190 : // leaving it to the user is because the value of the exponent that
191 : // should be used to decide if it is safe to switch back to DG should be
192 : // `persson_exponent+1` to prevent the code from rapidly switching back
193 : // and forth between DG and subcell. Rather than trying to enforce this
194 : // by documentation, the switching back to DG TCI gets passed in the
195 : // exponent it should use, and to keep the interface between the TCIs
196 : // consistent, we also pass the exponent in separately here.
197 : std::tuple<int, RdmpTciData> tci_result = db::mutate_apply<TciMutator>(
198 : make_not_null(&box), subcell_options.persson_exponent(), false);
199 :
200 : const int tci_decision = std::get<0>(tci_result);
201 : db::mutate<Tags::TciDecision>(
202 : [&tci_decision](const gsl::not_null<int*> tci_decision_ptr) {
203 : *tci_decision_ptr = tci_decision;
204 : },
205 : make_not_null(&box));
206 :
207 : cell_is_troubled |= (tci_decision != 0);
208 :
209 : // If either:
210 : //
211 : // 1. we are not allowed to do subcell in this block (handled by the
212 : // early return above)
213 : // 2. the element is at an outer boundary _and_ we aren't allowed to go
214 : // to subcell at an outer boundary.
215 : // 3. the cell is not troubled
216 : //
217 : // then we can remove the current neighbor data and update the RDMP TCI
218 : // data.
219 : if ((cell_has_external_boundary and
220 : not subcell_enabled_at_external_boundary) or
221 : not cell_is_troubled) {
222 : db::mutate<subcell::Tags::GhostDataForReconstruction<Dim>,
223 : subcell::Tags::DataForRdmpTci>(
224 : [&tci_result](const auto neighbor_data_ptr,
225 : const gsl::not_null<RdmpTciData*> rdmp_tci_data_ptr) {
226 : neighbor_data_ptr->clear();
227 : *rdmp_tci_data_ptr = std::move(std::get<1>(std::move(tci_result)));
228 : },
229 : make_not_null(&box));
230 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
231 : }
232 :
233 : db::mutate<variables_tag, ::Tags::HistoryEvolvedVariables<variables_tag>,
234 : Tags::ActiveGrid, Tags::DidRollback,
235 : subcell::Tags::GhostDataForReconstruction<Dim>>(
236 : [&dg_mesh, &element, &subcell_mesh](
237 : const auto active_vars_ptr, const auto active_history_ptr,
238 : const gsl::not_null<ActiveGrid*> active_grid_ptr,
239 : const gsl::not_null<bool*> did_rollback_ptr,
240 : const gsl::not_null<DirectionalIdMap<Dim, GhostData>*>
241 : ghost_data_ptr,
242 : const DirectionalIdMap<Dim, Mesh<Dim>>& meshes_for_ghost_data,
243 : const size_t ghost_zone_size,
244 : const DirectionalIdMap<Dim, std::optional<intrp::Irregular<Dim>>>&
245 : neighbor_dg_to_fd_interpolants) {
246 : ASSERT(active_history_ptr->size() > 0,
247 : "We cannot have an empty history when unwinding, that's just "
248 : "nutty. Did you call the action too early in the action "
249 : "list?");
250 : // Rollback u^{n+1}* to u^n (undoing the candidate solution).
251 : //
252 : // Note: strictly speaking, to be conservative this should project
253 : // uJ instead of u.
254 : *active_vars_ptr = fd::project(active_history_ptr->latest_value(),
255 : dg_mesh, subcell_mesh.extents());
256 :
257 : // Project the time stepper history to the subcells, excluding the
258 : // most recent inadmissible history.
259 : active_history_ptr->undo_latest();
260 : active_history_ptr->map_entries(
261 : [&dg_mesh, &subcell_mesh](const auto entry) {
262 : *entry = fd::project(*entry, dg_mesh, subcell_mesh.extents());
263 : });
264 : *active_grid_ptr = ActiveGrid::Subcell;
265 : *did_rollback_ptr = true;
266 : // Project the neighbor data we were sent for reconstruction since
267 : // the neighbor might have sent DG volume data instead of ghost data
268 : // in order to elide projections when they aren't necessary.
269 : for (const auto& [directional_element_id, mesh_for_ghost_data] :
270 : meshes_for_ghost_data) {
271 : evolution::dg::subcell::insert_or_update_neighbor_volume_data<
272 : false>(ghost_data_ptr,
273 : ghost_data_ptr->at(directional_element_id)
274 : .neighbor_ghost_data_for_reconstruction(),
275 : 0, directional_element_id, mesh_for_ghost_data, element,
276 : subcell_mesh, ghost_zone_size,
277 : neighbor_dg_to_fd_interpolants,
278 : typename Metavariables::SubcellOptions::GhostVariables::
279 : ghost_variables_tag_list{});
280 : }
281 :
282 : // Note: We do _not_ project the boundary history here because
283 : // that needs to be done at the lifting stage of the subcell
284 : // method, since we need to lift G+D instead of the ingredients
285 : // that go into G+D, which is what we would be projecting here.
286 : },
287 : make_not_null(&box),
288 : db::get<evolution::dg::subcell::Tags::MeshForGhostData<Dim>>(box),
289 : Metavariables::SubcellOptions::ghost_zone_size(box),
290 : db::get<
291 : evolution::dg::subcell::Tags::InterpolatorsFromNeighborDgToFd<Dim>>(
292 : box));
293 :
294 : if (UNLIKELY(db::get<::Tags::TimeStepId>(box).slab_number() < 0)) {
295 : // If we are doing self start, then we need to project the initial
296 : // guess to the subcells as well.
297 : //
298 : // Warning: this unfortunately needs to be kept in sync with the
299 : // self-start procedure.
300 : //
301 : // Note: if we switch to the subcells then we might have an
302 : // inconsistent
303 : // state between the primitive and conservative variables on the
304 : // subcells. The most correct thing is to re-compute the
305 : // primitive variables on the subcells, since projecting the
306 : // conservative variables is conservative.
307 : if constexpr (Metavariables::system::
308 : has_primitive_and_conservative_vars) {
309 : db::mutate<
310 : SelfStart::Tags::InitialValue<variables_tag>,
311 : SelfStart::Tags::InitialValue<
312 : typename Metavariables::system::primitive_variables_tag>>(
313 : [&dg_mesh, &subcell_mesh](const auto initial_vars_ptr,
314 : const auto initial_prim_vars_ptr) {
315 : // Note: for strict conservation, we need to project uJ
316 : // instead of just u.
317 : std::get<0>(*initial_vars_ptr) =
318 : fd::project(std::get<0>(*initial_vars_ptr), dg_mesh,
319 : subcell_mesh.extents());
320 : std::get<0>(*initial_prim_vars_ptr) =
321 : fd::project(std::get<0>(*initial_prim_vars_ptr), dg_mesh,
322 : subcell_mesh.extents());
323 : },
324 : make_not_null(&box));
325 : } else {
326 : db::mutate<SelfStart::Tags::InitialValue<variables_tag>>(
327 : [&dg_mesh, &subcell_mesh](const auto initial_vars_ptr) {
328 : // Note: for strict conservation, we need to project uJ
329 : // instead of just u.
330 : std::get<0>(*initial_vars_ptr) =
331 : fd::project(std::get<0>(*initial_vars_ptr), dg_mesh,
332 : subcell_mesh.extents());
333 : },
334 : make_not_null(&box));
335 : }
336 : }
337 :
338 : return {Parallel::AlgorithmExecution::Continue,
339 : tmpl::index_of<
340 : ActionList,
341 : ::Actions::Label<evolution::dg::subcell::Actions::Labels::
342 : BeginSubcellAfterDgRollback>>::value +
343 : 1};
344 : }
345 : };
346 : } // namespace evolution::dg::subcell::Actions
|
