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 <tuple>
10 : #include <unordered_map>
11 : #include <utility>
12 : #include <vector>
13 :
14 : #include "DataStructures/DataBox/DataBox.hpp"
15 : #include "DataStructures/DataBox/DataBoxTag.hpp"
16 : #include "DataStructures/DataBox/Prefixes.hpp"
17 : #include "DataStructures/Variables.hpp"
18 : #include "Domain/Creators/Tags/Domain.hpp"
19 : #include "Domain/Structure/ChildSize.hpp"
20 : #include "Domain/Structure/Direction.hpp"
21 : #include "Domain/Structure/Element.hpp"
22 : #include "Domain/Structure/Neighbors.hpp"
23 : #include "Domain/Structure/OrientationMap.hpp"
24 : #include "Domain/Structure/TrimMap.hpp"
25 : #include "Domain/Tags.hpp"
26 : #include "Domain/Tags/NeighborMesh.hpp"
27 : #include "Evolution/DiscontinuousGalerkin/InboxTags.hpp"
28 : #include "Evolution/DiscontinuousGalerkin/Initialization/QuadratureTag.hpp"
29 : #include "Evolution/DiscontinuousGalerkin/MortarData.hpp"
30 : #include "Evolution/DiscontinuousGalerkin/MortarDataHolder.hpp"
31 : #include "Evolution/DiscontinuousGalerkin/MortarInfo.hpp"
32 : #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
33 : #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp"
34 : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
35 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
36 : #include "NumericalAlgorithms/Spectral/SegmentSize.hpp"
37 : #include "Parallel/AlgorithmExecution.hpp"
38 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
39 : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
40 : #include "Time/BoundaryHistory.hpp"
41 : #include "Time/LtsMode.hpp"
42 : #include "Time/TimeStepId.hpp"
43 : #include "Utilities/ErrorHandling/Assert.hpp"
44 : #include "Utilities/Gsl.hpp"
45 : #include "Utilities/MakeArray.hpp"
46 : #include "Utilities/TMPL.hpp"
47 :
48 : /// \cond
49 : template <size_t Dim>
50 : class Domain;
51 : namespace Parallel {
52 : template <typename Metavariables>
53 : class GlobalCache;
54 : } // namespace Parallel
55 : namespace Spectral {
56 : enum class Quadrature : uint8_t;
57 : } // namespace Spectral
58 : namespace Tags {
59 : struct LtsMode;
60 : struct TimeStepId;
61 : } // namespace Tags
62 : namespace tuples {
63 : template <class... Tags>
64 : class TaggedTuple;
65 : } // namespace tuples
66 : /// \endcond
67 :
68 : namespace evolution::dg::Initialization {
69 : namespace detail {
70 : template <size_t Dim>
71 : ::dg::MortarMap<Dim, evolution::dg::MortarDataHolder<Dim>> empty_mortar_data(
72 : const Element<Dim>& element);
73 :
74 : template <size_t Dim>
75 : ::dg::MortarMap<Dim, MortarInfo<Dim>> mortar_infos(
76 : const Domain<Dim>& domain, const Element<Dim>& element,
77 : const Mesh<Dim>& volume_mesh,
78 : const ::dg::MortarMap<Dim, Mesh<Dim>>& neighbor_mesh, LtsMode lts_mode);
79 :
80 : template <size_t Dim>
81 : std::tuple<::dg::MortarMap<Dim, Mesh<Dim - 1>>,
82 : ::dg::MortarMap<Dim, TimeStepId>,
83 : DirectionMap<Dim, std::optional<Variables<tmpl::list<
84 : evolution::dg::Tags::MagnitudeOfNormal,
85 : evolution::dg::Tags::NormalCovector<Dim>>>>>>
86 : mortars_apply_impl(const Element<Dim>& element,
87 : const TimeStepId& next_temporal_id,
88 : const Mesh<Dim>& volume_mesh,
89 : const ::dg::MortarMap<Dim, Mesh<Dim>>& neighbor_mesh);
90 :
91 : template <size_t Dim>
92 : void h_refine_structure(
93 : gsl::not_null<::dg::MortarMap<Dim, evolution::dg::MortarDataHolder<Dim>>*>
94 : mortar_data,
95 : gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim - 1>>*> mortar_mesh,
96 : gsl::not_null<::dg::MortarMap<Dim, MortarInfo<Dim>>*> mortar_infos,
97 : gsl::not_null<::dg::MortarMap<Dim, TimeStepId>*> mortar_next_temporal_id,
98 : gsl::not_null<
99 : DirectionMap<Dim, std::optional<::Variables<tmpl::list<
100 : ::evolution::dg::Tags::MagnitudeOfNormal,
101 : ::evolution::dg::Tags::NormalCovector<Dim>>>>>*>
102 : normal_covector_and_magnitude,
103 : const Domain<Dim>& domain, const Mesh<Dim>& new_mesh,
104 : const Element<Dim>& new_element,
105 : const ::dg::MortarMap<Dim, Mesh<Dim>>& neighbor_mesh,
106 : const TimeStepId& current_temporal_id, LtsMode lts_mode);
107 : } // namespace detail
108 :
109 : /*!
110 : * \brief Initialize mortars between elements for exchanging boundary correction
111 : * terms.
112 : *
113 : * Uses:
114 : * - DataBox:
115 : * - `Tags::Element<Dim>`
116 : * - `Tags::Mesh<Dim>`
117 : * - `BoundaryScheme::receive_temporal_id`
118 : *
119 : * DataBox changes:
120 : * - Adds:
121 : * - `Tags::MortarData<Dim>`
122 : * - `Tags::MortarMesh<Dim>`
123 : * - `Tags::MortarInfo<Dim>`
124 : * - `Tags::MortarNextTemporalId<Dim>`
125 : * - `evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>`
126 : * - Removes: nothing
127 : * - Modifies: nothing
128 : */
129 : template <size_t Dim>
130 1 : struct Mortars {
131 : public:
132 0 : using const_global_cache_tags = tmpl::list<domain::Tags::Domain<Dim>>;
133 0 : using simple_tags_from_options = tmpl::list<>;
134 :
135 0 : using simple_tags =
136 : tmpl::list<Tags::MortarData<Dim>, Tags::MortarMesh<Dim>,
137 : Tags::MortarInfo<Dim>, Tags::MortarNextTemporalId<Dim>,
138 : evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>,
139 : Tags::MortarDataHistory<Dim>>;
140 0 : using compute_tags = tmpl::list<>;
141 :
142 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
143 : typename ArrayIndex, typename ActionList,
144 : typename ParallelComponent>
145 0 : static Parallel::iterable_action_return_t apply(
146 : db::DataBox<DbTagsList>& box,
147 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
148 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
149 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
150 : const ParallelComponent* const /*meta*/) {
151 : const auto& domain = db::get<domain::Tags::Domain<Dim>>(box);
152 : const auto& element = db::get<::domain::Tags::Element<Dim>>(box);
153 : const auto& volume_mesh = db::get<domain::Tags::Mesh<Dim>>(box);
154 : const auto& neighbor_mesh = db::get<domain::Tags::NeighborMesh<Dim>>(box);
155 : const auto lts_mode = db::get<::Tags::LtsMode>(box);
156 : auto mortar_data = detail::empty_mortar_data(element);
157 : auto mortar_infos = detail::mortar_infos(domain, element, volume_mesh,
158 : neighbor_mesh, lts_mode);
159 : auto [mortar_meshes, mortar_next_temporal_ids, normal_covector_quantities] =
160 : detail::mortars_apply_impl(
161 : element, db::get<::Tags::Next<::Tags::TimeStepId>>(box),
162 : db::get<::domain::Tags::Mesh<Dim>>(box),
163 : db::get<::domain::Tags::NeighborMesh<Dim>>(box));
164 : typename Tags::MortarDataHistory<Dim>::type boundary_data_history{};
165 : for (const auto& mortar_id_and_data : mortar_data) {
166 : if (mortar_infos.at(mortar_id_and_data.first).time_stepping_policy() ==
167 : LtsMode::Conservative) {
168 : // default initialize data
169 : boundary_data_history[mortar_id_and_data.first];
170 : }
171 : }
172 : ::Initialization::mutate_assign<simple_tags>(
173 : make_not_null(&box), std::move(mortar_data), std::move(mortar_meshes),
174 : std::move(mortar_infos), std::move(mortar_next_temporal_ids),
175 : std::move(normal_covector_quantities),
176 : std::move(boundary_data_history));
177 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
178 : }
179 : };
180 :
181 : /// \brief Initialize/update items related to mortars after an AMR change
182 : ///
183 : /// Mutates:
184 : /// - Tags::MortarData<dim>
185 : /// - Tags::MortarMesh<dim>
186 : /// - Tags::MortarInfo<dim>
187 : /// - Tags::MortarNextTemporalId<dim>
188 : /// - evolution::dg::Tags::NormalCovectorAndMagnitude<dim>
189 : /// - Tags::MortarDataHistory<dim>>
190 : ///
191 : /// For p-refined interfaces:
192 : /// - Regenerates MortarData and MortarInfo (should have no effect)
193 : /// - Sets the NormalCovectorAndMagnitude to std::nullopt if the face mesh
194 : /// changed
195 : /// - Projects the local geometric data (but not the data on the mortar-mesh)
196 : /// in the MortarDataHistory, if present
197 : /// - Does nothing to MortarMesh and MortarNextTemporalId
198 : ///
199 : /// For h-refined interfaces:
200 : /// - Regenerates MortarData and MortarInfo
201 : /// - Sets the NormalCovectorAndMagnitude to std::nullopt
202 : /// - Calculates MortarMesh
203 : /// - Sets MortarNextTemporalId to the current temporal id
204 : /// - For local time-stepping:
205 : /// - Removes MortarDataHistory data corresponding to split or joined
206 : /// elements
207 : /// - Projects MortarDataHistory data corresponding to non-h-refined
208 : /// elements onto refined mortars (both geometric and mortar-mesh data)
209 : /// - Creates empty histories for new mortars between two h-refined
210 : /// elements
211 : template <size_t Dim>
212 1 : struct ProjectMortars : tt::ConformsTo<amr::protocols::Projector> {
213 : private:
214 0 : using magnitude_and_normal_type =
215 : ::Variables<tmpl::list<::evolution::dg::Tags::MagnitudeOfNormal,
216 : ::evolution::dg::Tags::NormalCovector<Dim>>>;
217 :
218 : public:
219 0 : using mortar_data_history_tag = Tags::MortarDataHistory<Dim>;
220 0 : using mortar_data_history_type = typename mortar_data_history_tag::type;
221 :
222 0 : using return_tags =
223 : tmpl::list<Tags::MortarData<Dim>, Tags::MortarMesh<Dim>,
224 : Tags::MortarInfo<Dim>, Tags::MortarNextTemporalId<Dim>,
225 : evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>,
226 : Tags::MortarDataHistory<Dim>>;
227 0 : using argument_tags =
228 : tmpl::list<domain::Tags::Domain<Dim>, domain::Tags::Mesh<Dim>,
229 : domain::Tags::Element<Dim>, domain::Tags::NeighborMesh<Dim>,
230 : ::Tags::TimeStepId, ::Tags::LtsMode>;
231 :
232 0 : static void apply(
233 : gsl::not_null<::dg::MortarMap<Dim, evolution::dg::MortarDataHolder<Dim>>*>
234 : mortar_data,
235 : gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim - 1>>*> mortar_mesh,
236 : gsl::not_null<::dg::MortarMap<Dim, MortarInfo<Dim>>*> mortar_infos,
237 : gsl::not_null<::dg::MortarMap<Dim, TimeStepId>*> mortar_next_temporal_id,
238 : gsl::not_null<
239 : DirectionMap<Dim, std::optional<magnitude_and_normal_type>>*>
240 : normal_covector_and_magnitude,
241 : gsl::not_null<mortar_data_history_type*> mortar_data_history,
242 : const Domain<Dim>& domain, const Mesh<Dim>& new_mesh,
243 : const Element<Dim>& new_element,
244 : const ::dg::MortarMap<Dim, Mesh<Dim>>& neighbor_mesh,
245 : const TimeStepId& current_temporal_id, LtsMode lts_mode,
246 : const std::pair<Mesh<Dim>, Element<Dim>>& old_mesh_and_element);
247 :
248 : template <typename... ParentTags>
249 0 : static void apply(
250 : const gsl::not_null<
251 : ::dg::MortarMap<Dim, evolution::dg::MortarDataHolder<Dim>>*>
252 : mortar_data,
253 : const gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim - 1>>*> mortar_mesh,
254 : const gsl::not_null<::dg::MortarMap<Dim, MortarInfo<Dim>>*> mortar_infos,
255 : const gsl::not_null<::dg::MortarMap<Dim, TimeStepId>*>
256 : mortar_next_temporal_id,
257 : const gsl::not_null<
258 : DirectionMap<Dim, std::optional<magnitude_and_normal_type>>*>
259 : normal_covector_and_magnitude,
260 : const gsl::not_null<mortar_data_history_type*> mortar_data_history,
261 : const Domain<Dim>& domain, const Mesh<Dim>& new_mesh,
262 : const Element<Dim>& new_element,
263 : const ::dg::MortarMap<Dim, Mesh<Dim>>& neighbor_mesh,
264 : const TimeStepId& /*possibly_unset*/, const LtsMode lts_mode,
265 : const tuples::TaggedTuple<ParentTags...>& parent_items) {
266 : detail::h_refine_structure(
267 : mortar_data, mortar_mesh, mortar_infos, mortar_next_temporal_id,
268 : normal_covector_and_magnitude, domain, new_mesh, new_element,
269 : neighbor_mesh, get<::Tags::TimeStepId>(parent_items), lts_mode);
270 :
271 : const auto& old_element = get<domain::Tags::Element<Dim>>(parent_items);
272 : const auto& old_histories = get<mortar_data_history_tag>(parent_items);
273 : for (const auto& [direction, neighbors] : new_element.neighbors()) {
274 : for (const auto& neighbor : neighbors) {
275 : const DirectionalId<Dim> mortar_id{direction, neighbor};
276 : if (mortar_infos->at(mortar_id).time_stepping_policy() !=
277 : LtsMode::Conservative) {
278 : continue;
279 : }
280 : if (const auto old_history = old_histories.find(mortar_id);
281 : old_history != old_histories.end()) {
282 : // The neighbor did not h-refine, so we have to project
283 : // its mortar data from our parent.
284 : auto& new_history =
285 : mortar_data_history->emplace(mortar_id, old_history->second)
286 : .first->second;
287 : new_history.local().clear();
288 : auto remote_history = new_history.remote();
289 : const auto& new_mortar_mesh = mortar_mesh->at(mortar_id);
290 : const auto& orientation = neighbors.orientation(neighbor);
291 : const auto new_mortar_size = domain::child_size(
292 : ::dg::mortar_segments(new_element.id(), neighbor,
293 : direction.dimension(), orientation),
294 : ::dg::mortar_segments(old_element.id(), neighbor,
295 : direction.dimension(), orientation));
296 : const auto project_mortar_data =
297 : [&new_mortar_mesh, &new_mortar_size](
298 : const TimeStepId& /* id */,
299 : const gsl::not_null<::evolution::dg::MortarData<Dim>*> data) {
300 : const auto& old_mortar_mesh = data->mortar_mesh.value();
301 : DataVector& vars = data->mortar_data.value();
302 : vars = Spectral::project(
303 : vars, old_mortar_mesh, new_mortar_mesh,
304 : make_array<Dim - 1>(Spectral::SegmentSize::Full),
305 : new_mortar_size);
306 : data->mortar_mesh = new_mortar_mesh;
307 : return true;
308 : };
309 : remote_history.for_each(project_mortar_data);
310 : } else {
311 : // Neither this element nor the neighbor existed before
312 : // refinement.
313 : mortar_data_history->emplace(
314 : mortar_id, typename mortar_data_history_type::mapped_type{});
315 : }
316 : }
317 : }
318 : }
319 :
320 : template <typename... ChildTags>
321 0 : static void apply(
322 : const gsl::not_null<
323 : ::dg::MortarMap<Dim, evolution::dg::MortarDataHolder<Dim>>*>
324 : mortar_data,
325 : const gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim - 1>>*> mortar_mesh,
326 : const gsl::not_null<::dg::MortarMap<Dim, MortarInfo<Dim>>*> mortar_infos,
327 : const gsl::not_null<::dg::MortarMap<Dim, TimeStepId>*>
328 : mortar_next_temporal_id,
329 : const gsl::not_null<
330 : DirectionMap<Dim, std::optional<magnitude_and_normal_type>>*>
331 : normal_covector_and_magnitude,
332 : const gsl::not_null<mortar_data_history_type*> mortar_data_history,
333 : const Domain<Dim>& domain, const Mesh<Dim>& new_mesh,
334 : const Element<Dim>& new_element,
335 : const ::dg::MortarMap<Dim, Mesh<Dim>>& neighbor_mesh,
336 : const TimeStepId& /*possibly_unset*/, const LtsMode lts_mode,
337 : const std::unordered_map<
338 : ElementId<Dim>, tuples::TaggedTuple<ChildTags...>>& children_items) {
339 : detail::h_refine_structure(
340 : mortar_data, mortar_mesh, mortar_infos, mortar_next_temporal_id,
341 : normal_covector_and_magnitude, domain, new_mesh, new_element,
342 : neighbor_mesh, get<::Tags::TimeStepId>(children_items.begin()->second),
343 : lts_mode);
344 :
345 : for (const auto& [direction, neighbors] : new_element.neighbors()) {
346 : for (const auto& neighbor : neighbors) {
347 : const DirectionalId<Dim> mortar_id{direction, neighbor};
348 : if (mortar_infos->at(mortar_id).time_stepping_policy() !=
349 : LtsMode::Conservative) {
350 : continue;
351 : }
352 : std::optional<typename mortar_data_history_type::mapped_type>
353 : new_history{};
354 : for (const auto& [child, child_items] : children_items) {
355 : const auto& old_histories =
356 : get<mortar_data_history_tag>(child_items);
357 : if (const auto old_history = old_histories.find(mortar_id);
358 : old_history != old_histories.end()) {
359 : // The neighbor did not h-refine, so we have to project
360 : // its mortar data from our children.
361 : const auto& new_mortar_mesh = mortar_mesh->at(mortar_id);
362 : const auto& orientation = neighbors.orientation(neighbor);
363 : const auto old_mortar_size = domain::child_size(
364 : ::dg::mortar_segments(child, neighbor, direction.dimension(),
365 : orientation),
366 : ::dg::mortar_segments(new_element.id(), neighbor,
367 : direction.dimension(), orientation));
368 : if (not new_history.has_value()) {
369 : new_history.emplace(old_history->second);
370 : new_history->local().clear();
371 : auto remote_history = new_history->remote();
372 : const auto project_mortar_data =
373 : [&new_mortar_mesh, &old_mortar_size](
374 : const TimeStepId& /* id */,
375 : const gsl::not_null<::evolution::dg::MortarData<Dim>*>
376 : data) {
377 : const auto& old_mortar_mesh = data->mortar_mesh.value();
378 : DataVector& vars = data->mortar_data.value();
379 : vars = Spectral::project(
380 : vars, old_mortar_mesh, new_mortar_mesh, old_mortar_size,
381 : make_array<Dim - 1>(Spectral::SegmentSize::Full));
382 : data->mortar_mesh = new_mortar_mesh;
383 : return true;
384 : };
385 : remote_history.for_each(project_mortar_data);
386 : } else {
387 : auto remote_history = new_history->remote();
388 : const auto old_remote_history = old_history->second.remote();
389 : const auto project_mortar_data =
390 : [&new_mortar_mesh, &old_mortar_size, &old_remote_history](
391 : const TimeStepId& id,
392 : const gsl::not_null<::evolution::dg::MortarData<Dim>*>
393 : data) {
394 : const auto& old_data = old_remote_history.data(id);
395 : const auto& old_mortar_mesh =
396 : old_data.mortar_mesh.value();
397 : data->mortar_data.value() += Spectral::project(
398 : old_data.mortar_data.value(), old_mortar_mesh,
399 : new_mortar_mesh, old_mortar_size,
400 : make_array<Dim - 1>(Spectral::SegmentSize::Full));
401 : return true;
402 : };
403 : remote_history.for_each(project_mortar_data);
404 : }
405 : }
406 : }
407 :
408 : if (new_history.has_value()) {
409 : mortar_data_history->emplace(mortar_id, std::move(*new_history));
410 : } else {
411 : // Neither this element nor the neighbor existed before
412 : // refinement.
413 : mortar_data_history->emplace(
414 : mortar_id, typename mortar_data_history_type::mapped_type{});
415 : }
416 : }
417 : }
418 : }
419 : };
420 : } // namespace evolution::dg::Initialization
|
