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/ApplyMatrices.hpp"
15 : #include "DataStructures/DataBox/DataBox.hpp"
16 : #include "DataStructures/DataBox/DataBoxTag.hpp"
17 : #include "DataStructures/DataBox/Prefixes.hpp"
18 : #include "DataStructures/Variables.hpp"
19 : #include "Domain/Creators/Tags/Domain.hpp"
20 : #include "Domain/Structure/ChildSize.hpp"
21 : #include "Domain/Structure/Direction.hpp"
22 : #include "Domain/Structure/Element.hpp"
23 : #include "Domain/Structure/Neighbors.hpp"
24 : #include "Domain/Structure/OrientationMap.hpp"
25 : #include "Domain/Structure/TrimMap.hpp"
26 : #include "Domain/Tags.hpp"
27 : #include "Domain/Tags/NeighborMesh.hpp"
28 : #include "Evolution/DiscontinuousGalerkin/InboxTags.hpp"
29 : #include "Evolution/DiscontinuousGalerkin/Initialization/QuadratureTag.hpp"
30 : #include "Evolution/DiscontinuousGalerkin/MortarData.hpp"
31 : #include "Evolution/DiscontinuousGalerkin/MortarDataHolder.hpp"
32 : #include "Evolution/DiscontinuousGalerkin/MortarInfo.hpp"
33 : #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
34 : #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp"
35 : #include "Evolution/DiscontinuousGalerkin/TimeSteppingPolicy.hpp"
36 : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
37 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
38 : #include "NumericalAlgorithms/Spectral/SegmentSize.hpp"
39 : #include "Parallel/AlgorithmExecution.hpp"
40 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
41 : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
42 : #include "Time/BoundaryHistory.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 TimeStepId;
61 : } // namespace Tags
62 : namespace tuples {
63 : template <class... Tags>
64 : class TaggedTuple;
65 : } // namespace tuples
66 : /// \endcond
67 :
68 1 : 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,
79 : bool local_time_stepping);
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, bool local_time_stepping);
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 : auto mortar_data = detail::empty_mortar_data(element);
157 : auto mortar_infos =
158 : detail::mortar_infos(domain, element, volume_mesh, neighbor_mesh,
159 : Metavariables::local_time_stepping);
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, bool LocalTimeStepping>
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>;
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,
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*/,
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),
271 : LocalTimeStepping);
272 :
273 : const auto& old_element = get<domain::Tags::Element<Dim>>(parent_items);
274 : const auto& old_histories = get<mortar_data_history_tag>(parent_items);
275 : for (const auto& [direction, neighbors] : new_element.neighbors()) {
276 : for (const auto& neighbor : neighbors) {
277 : const DirectionalId<Dim> mortar_id{direction, neighbor};
278 : if (mortar_infos->at(mortar_id).time_stepping_policy() !=
279 : TimeSteppingPolicy::Conservative) {
280 : continue;
281 : }
282 : if (const auto old_history = old_histories.find(mortar_id);
283 : old_history != old_histories.end()) {
284 : // The neighbor did not h-refine, so we have to project
285 : // its mortar data from our parent.
286 : auto& new_history =
287 : mortar_data_history->emplace(mortar_id, old_history->second)
288 : .first->second;
289 : new_history.local().clear();
290 : auto remote_history = new_history.remote();
291 : const auto& new_mortar_mesh = mortar_mesh->at(mortar_id);
292 : const auto& orientation = neighbors.orientation(neighbor);
293 : const auto new_mortar_size = domain::child_size(
294 : ::dg::mortar_segments(new_element.id(), neighbor,
295 : direction.dimension(), orientation),
296 : ::dg::mortar_segments(old_element.id(), neighbor,
297 : direction.dimension(), orientation));
298 : const auto project_mortar_data =
299 : [&new_mortar_mesh, &new_mortar_size](
300 : const TimeStepId& /* id */,
301 : const gsl::not_null<::evolution::dg::MortarData<Dim>*> data) {
302 : const auto& old_mortar_mesh = data->mortar_mesh.value();
303 : const auto mortar_projection_matrices =
304 : Spectral::projection_matrices(
305 : old_mortar_mesh, new_mortar_mesh,
306 : make_array<Dim - 1>(Spectral::SegmentSize::Full),
307 : new_mortar_size);
308 : DataVector& vars = data->mortar_data.value();
309 : vars = apply_matrices(mortar_projection_matrices, vars,
310 : old_mortar_mesh.extents());
311 : data->mortar_mesh = new_mortar_mesh;
312 : return true;
313 : };
314 : remote_history.for_each(project_mortar_data);
315 : } else {
316 : // Neither this element nor the neighbor existed before
317 : // refinement.
318 : mortar_data_history->emplace(
319 : mortar_id, typename mortar_data_history_type::mapped_type{});
320 : }
321 : }
322 : }
323 : }
324 :
325 : template <typename... ChildTags>
326 0 : static void apply(
327 : const gsl::not_null<
328 : ::dg::MortarMap<Dim, evolution::dg::MortarDataHolder<Dim>>*>
329 : mortar_data,
330 : const gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim - 1>>*> mortar_mesh,
331 : const gsl::not_null<::dg::MortarMap<Dim, MortarInfo<Dim>>*> mortar_infos,
332 : const gsl::not_null<::dg::MortarMap<Dim, TimeStepId>*>
333 : mortar_next_temporal_id,
334 : const gsl::not_null<
335 : DirectionMap<Dim, std::optional<magnitude_and_normal_type>>*>
336 : normal_covector_and_magnitude,
337 : const gsl::not_null<mortar_data_history_type*> mortar_data_history,
338 : const Domain<Dim>& domain, const Mesh<Dim>& new_mesh,
339 : const Element<Dim>& new_element,
340 : const ::dg::MortarMap<Dim, Mesh<Dim>>& neighbor_mesh,
341 : const TimeStepId& /*possibly_unset*/,
342 : const std::unordered_map<
343 : ElementId<Dim>, tuples::TaggedTuple<ChildTags...>>& children_items) {
344 : detail::h_refine_structure(
345 : mortar_data, mortar_mesh, mortar_infos, mortar_next_temporal_id,
346 : normal_covector_and_magnitude, domain, new_mesh, new_element,
347 : neighbor_mesh, get<::Tags::TimeStepId>(children_items.begin()->second),
348 : LocalTimeStepping);
349 :
350 : for (const auto& [direction, neighbors] : new_element.neighbors()) {
351 : for (const auto& neighbor : neighbors) {
352 : const DirectionalId<Dim> mortar_id{direction, neighbor};
353 : if (mortar_infos->at(mortar_id).time_stepping_policy() !=
354 : TimeSteppingPolicy::Conservative) {
355 : continue;
356 : }
357 : std::optional<typename mortar_data_history_type::mapped_type>
358 : new_history{};
359 : for (const auto& [child, child_items] : children_items) {
360 : const auto& old_histories =
361 : get<mortar_data_history_tag>(child_items);
362 : if (const auto old_history = old_histories.find(mortar_id);
363 : old_history != old_histories.end()) {
364 : // The neighbor did not h-refine, so we have to project
365 : // its mortar data from our children.
366 : const auto& new_mortar_mesh = mortar_mesh->at(mortar_id);
367 : const auto& orientation = neighbors.orientation(neighbor);
368 : const auto old_mortar_size = domain::child_size(
369 : ::dg::mortar_segments(child, neighbor, direction.dimension(),
370 : orientation),
371 : ::dg::mortar_segments(new_element.id(), neighbor,
372 : direction.dimension(), orientation));
373 : if (not new_history.has_value()) {
374 : new_history.emplace(old_history->second);
375 : new_history->local().clear();
376 : auto remote_history = new_history->remote();
377 : const auto project_mortar_data =
378 : [&new_mortar_mesh, &old_mortar_size](
379 : const TimeStepId& /* id */,
380 : const gsl::not_null<::evolution::dg::MortarData<Dim>*>
381 : data) {
382 : const auto& old_mortar_mesh = data->mortar_mesh.value();
383 : const auto mortar_projection_matrices =
384 : Spectral::projection_matrices(
385 : old_mortar_mesh, new_mortar_mesh, old_mortar_size,
386 : make_array<Dim - 1>(Spectral::SegmentSize::Full));
387 : DataVector& vars = data->mortar_data.value();
388 : vars = apply_matrices(mortar_projection_matrices, vars,
389 : old_mortar_mesh.extents());
390 : data->mortar_mesh = new_mortar_mesh;
391 : return true;
392 : };
393 : remote_history.for_each(project_mortar_data);
394 : } else {
395 : auto remote_history = new_history->remote();
396 : const auto old_remote_history = old_history->second.remote();
397 : const auto project_mortar_data =
398 : [&new_mortar_mesh, &old_mortar_size, &old_remote_history](
399 : const TimeStepId& id,
400 : const gsl::not_null<::evolution::dg::MortarData<Dim>*>
401 : data) {
402 : const auto& old_data = old_remote_history.data(id);
403 : const auto& old_mortar_mesh =
404 : old_data.mortar_mesh.value();
405 : const auto mortar_projection_matrices =
406 : Spectral::projection_matrices(
407 : old_mortar_mesh, new_mortar_mesh, old_mortar_size,
408 : make_array<Dim - 1>(Spectral::SegmentSize::Full));
409 : data->mortar_data.value() +=
410 : apply_matrices(mortar_projection_matrices,
411 : old_data.mortar_data.value(),
412 : old_mortar_mesh.extents());
413 : return true;
414 : };
415 : remote_history.for_each(project_mortar_data);
416 : }
417 : }
418 : }
419 :
420 : if (new_history.has_value()) {
421 : mortar_data_history->emplace(mortar_id, std::move(*new_history));
422 : } else {
423 : // Neither this element nor the neighbor existed before
424 : // refinement.
425 : mortar_data_history->emplace(
426 : mortar_id, typename mortar_data_history_type::mapped_type{});
427 : }
428 : }
429 : }
430 : }
431 : };
432 : } // namespace evolution::dg::Initialization
|
