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