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/PrefixHelpers.hpp"
18 : #include "DataStructures/DataBox/Prefixes.hpp"
19 : #include "DataStructures/Variables.hpp"
20 : #include "Domain/Creators/Tags/Domain.hpp"
21 : #include "Domain/Structure/ChildSize.hpp"
22 : #include "Domain/Structure/Direction.hpp"
23 : #include "Domain/Structure/Element.hpp"
24 : #include "Domain/Structure/Neighbors.hpp"
25 : #include "Domain/Structure/OrientationMap.hpp"
26 : #include "Domain/Structure/TrimMap.hpp"
27 : #include "Domain/Tags.hpp"
28 : #include "Domain/Tags/NeighborMesh.hpp"
29 : #include "Evolution/DiscontinuousGalerkin/InboxTags.hpp"
30 : #include "Evolution/DiscontinuousGalerkin/Initialization/QuadratureTag.hpp"
31 : #include "Evolution/DiscontinuousGalerkin/MortarData.hpp"
32 : #include "Evolution/DiscontinuousGalerkin/MortarDataHolder.hpp"
33 : #include "Evolution/DiscontinuousGalerkin/MortarInfo.hpp"
34 : #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
35 : #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp"
36 : #include "Evolution/DiscontinuousGalerkin/TimeSteppingPolicy.hpp"
37 : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
38 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
39 : #include "NumericalAlgorithms/Spectral/SegmentSize.hpp"
40 : #include "Parallel/AlgorithmExecution.hpp"
41 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
42 : #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
43 : #include "Time/BoundaryHistory.hpp"
44 : #include "Time/TimeStepId.hpp"
45 : #include "Utilities/ErrorHandling/Assert.hpp"
46 : #include "Utilities/Gsl.hpp"
47 : #include "Utilities/MakeArray.hpp"
48 : #include "Utilities/TMPL.hpp"
49 :
50 : /// \cond
51 : template <size_t Dim>
52 : class Domain;
53 : namespace Parallel {
54 : template <typename Metavariables>
55 : class GlobalCache;
56 : } // namespace Parallel
57 : namespace Spectral {
58 : enum class Quadrature : uint8_t;
59 : } // namespace Spectral
60 : namespace Tags {
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,
80 : bool local_time_stepping);
81 :
82 : template <size_t Dim>
83 : std::tuple<::dg::MortarMap<Dim, Mesh<Dim - 1>>,
84 : ::dg::MortarMap<Dim, TimeStepId>,
85 : DirectionMap<Dim, std::optional<Variables<tmpl::list<
86 : evolution::dg::Tags::MagnitudeOfNormal,
87 : evolution::dg::Tags::NormalCovector<Dim>>>>>>
88 : mortars_apply_impl(const Element<Dim>& element,
89 : const TimeStepId& next_temporal_id,
90 : const Mesh<Dim>& volume_mesh,
91 : const ::dg::MortarMap<Dim, Mesh<Dim>>& neighbor_mesh);
92 :
93 : template <size_t Dim>
94 : void h_refine_structure(
95 : gsl::not_null<::dg::MortarMap<Dim, evolution::dg::MortarDataHolder<Dim>>*>
96 : mortar_data,
97 : gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim - 1>>*> mortar_mesh,
98 : gsl::not_null<::dg::MortarMap<Dim, MortarInfo<Dim>>*> mortar_infos,
99 : gsl::not_null<::dg::MortarMap<Dim, TimeStepId>*> mortar_next_temporal_id,
100 : gsl::not_null<
101 : DirectionMap<Dim, std::optional<::Variables<tmpl::list<
102 : ::evolution::dg::Tags::MagnitudeOfNormal,
103 : ::evolution::dg::Tags::NormalCovector<Dim>>>>>*>
104 : normal_covector_and_magnitude,
105 : const Domain<Dim>& domain, const Mesh<Dim>& new_mesh,
106 : const Element<Dim>& new_element,
107 : const ::dg::MortarMap<Dim, Mesh<Dim>>& neighbor_mesh,
108 : const TimeStepId& current_temporal_id, bool local_time_stepping);
109 : } // namespace detail
110 :
111 : /*!
112 : * \brief Initialize mortars between elements for exchanging boundary correction
113 : * terms.
114 : *
115 : * Uses:
116 : * - DataBox:
117 : * - `Tags::Element<Dim>`
118 : * - `Tags::Mesh<Dim>`
119 : * - `BoundaryScheme::receive_temporal_id`
120 : *
121 : * DataBox changes:
122 : * - Adds:
123 : * - `Tags::MortarData<Dim>`
124 : * - `Tags::MortarMesh<Dim>`
125 : * - `Tags::MortarInfo<Dim>`
126 : * - `Tags::MortarNextTemporalId<Dim>`
127 : * - `evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>`
128 : * - Removes: nothing
129 : * - Modifies: nothing
130 : */
131 : template <size_t Dim, typename System>
132 1 : struct Mortars {
133 : public:
134 0 : using const_global_cache_tags = tmpl::list<domain::Tags::Domain<Dim>>;
135 0 : using simple_tags_from_options = tmpl::list<>;
136 :
137 0 : using simple_tags = tmpl::list<
138 : Tags::MortarData<Dim>, Tags::MortarMesh<Dim>, Tags::MortarInfo<Dim>,
139 : Tags::MortarNextTemporalId<Dim>,
140 : evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>,
141 : Tags::MortarDataHistory<
142 : Dim, typename db::add_tag_prefix<
143 : ::Tags::dt, typename System::variables_tag>::type>>;
144 0 : using compute_tags = tmpl::list<>;
145 :
146 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
147 : typename ArrayIndex, typename ActionList,
148 : typename ParallelComponent>
149 0 : static Parallel::iterable_action_return_t apply(
150 : db::DataBox<DbTagsList>& box,
151 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
152 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
153 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
154 : const ParallelComponent* const /*meta*/) {
155 : const auto& domain = db::get<domain::Tags::Domain<Dim>>(box);
156 : const auto& element = db::get<::domain::Tags::Element<Dim>>(box);
157 : const auto& volume_mesh = db::get<domain::Tags::Mesh<Dim>>(box);
158 : const auto& neighbor_mesh = db::get<domain::Tags::NeighborMesh<Dim>>(box);
159 : auto mortar_data = detail::empty_mortar_data(element);
160 : auto mortar_infos =
161 : detail::mortar_infos(domain, element, volume_mesh, neighbor_mesh,
162 : Metavariables::local_time_stepping);
163 : auto [mortar_meshes, mortar_next_temporal_ids, normal_covector_quantities] =
164 : detail::mortars_apply_impl(
165 : element, db::get<::Tags::Next<::Tags::TimeStepId>>(box),
166 : db::get<::domain::Tags::Mesh<Dim>>(box),
167 : db::get<::domain::Tags::NeighborMesh<Dim>>(box));
168 : typename Tags::MortarDataHistory<
169 : Dim, typename db::add_tag_prefix<
170 : ::Tags::dt, typename System::variables_tag>::type>::type
171 : boundary_data_history{};
172 : for (const auto& mortar_id_and_data : mortar_data) {
173 : if (mortar_infos.at(mortar_id_and_data.first).time_stepping_policy() ==
174 : TimeSteppingPolicy::Conservative) {
175 : // default initialize data
176 : boundary_data_history[mortar_id_and_data.first];
177 : }
178 : }
179 : ::Initialization::mutate_assign<simple_tags>(
180 : make_not_null(&box), std::move(mortar_data), std::move(mortar_meshes),
181 : std::move(mortar_infos), std::move(mortar_next_temporal_ids),
182 : std::move(normal_covector_quantities),
183 : std::move(boundary_data_history));
184 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
185 : }
186 : };
187 :
188 : /// \brief Initialize/update items related to mortars after an AMR change
189 : ///
190 : /// Mutates:
191 : /// - Tags::MortarData<dim>
192 : /// - Tags::MortarMesh<dim>
193 : /// - Tags::MortarInfo<dim>
194 : /// - Tags::MortarNextTemporalId<dim>
195 : /// - evolution::dg::Tags::NormalCovectorAndMagnitude<dim>
196 : /// - Tags::MortarDataHistory<dim, typename dt_variables_tag::type>>
197 : ///
198 : /// For p-refined interfaces:
199 : /// - Regenerates MortarData and MortarInfo (should have no effect)
200 : /// - Sets the NormalCovectorAndMagnitude to std::nullopt if the face mesh
201 : /// changed
202 : /// - Projects the local geometric data (but not the data on the mortar-mesh)
203 : /// in the MortarDataHistory, if present
204 : /// - Does nothing to MortarMesh and MortarNextTemporalId
205 : ///
206 : /// For h-refined interfaces:
207 : /// - Regenerates MortarData and MortarInfo
208 : /// - Sets the NormalCovectorAndMagnitude to std::nullopt
209 : /// - Calculates MortarMesh
210 : /// - Sets MortarNextTemporalId to the current temporal id
211 : /// - For local time-stepping:
212 : /// - Removes MortarDataHistory data corresponding to split or joined
213 : /// elements
214 : /// - Projects MortarDataHistory data corresponding to non-h-refined
215 : /// elements onto refined mortars (both geometric and mortar-mesh data)
216 : /// - Creates empty histories for new mortars between two h-refined
217 : /// elements
218 : template <typename Metavariables>
219 1 : struct ProjectMortars : tt::ConformsTo<amr::protocols::Projector> {
220 : private:
221 0 : using magnitude_and_normal_type = ::Variables<tmpl::list<
222 : ::evolution::dg::Tags::MagnitudeOfNormal,
223 : ::evolution::dg::Tags::NormalCovector<Metavariables::volume_dim>>>;
224 :
225 : public:
226 0 : static constexpr size_t dim = Metavariables::volume_dim;
227 0 : using dt_variables_tag = typename db::add_tag_prefix<
228 : ::Tags::dt, typename Metavariables::system::variables_tag>;
229 0 : using mortar_data_history_tag =
230 : Tags::MortarDataHistory<dim, typename dt_variables_tag::type>;
231 0 : using mortar_data_history_type = typename mortar_data_history_tag::type;
232 :
233 0 : using return_tags =
234 : tmpl::list<Tags::MortarData<dim>, Tags::MortarMesh<dim>,
235 : Tags::MortarInfo<dim>, Tags::MortarNextTemporalId<dim>,
236 : evolution::dg::Tags::NormalCovectorAndMagnitude<dim>,
237 : Tags::MortarDataHistory<dim, typename dt_variables_tag::type>>;
238 0 : using argument_tags =
239 : tmpl::list<domain::Tags::Domain<dim>, domain::Tags::Mesh<dim>,
240 : domain::Tags::Element<dim>, domain::Tags::NeighborMesh<dim>,
241 : ::Tags::TimeStepId>;
242 :
243 0 : static void apply(
244 : const gsl::not_null<
245 : ::dg::MortarMap<dim, evolution::dg::MortarDataHolder<dim>>*>
246 : mortar_data,
247 : const gsl::not_null<::dg::MortarMap<dim, Mesh<dim - 1>>*> mortar_mesh,
248 : const gsl::not_null<::dg::MortarMap<dim, MortarInfo<dim>>*> mortar_infos,
249 : const gsl::not_null<::dg::MortarMap<dim, TimeStepId>*>
250 : mortar_next_temporal_id,
251 : const gsl::not_null<
252 : DirectionMap<dim, std::optional<magnitude_and_normal_type>>*>
253 : normal_covector_and_magnitude,
254 : const gsl::not_null<mortar_data_history_type*> mortar_data_history,
255 : const Domain<dim>& domain, const Mesh<dim>& new_mesh,
256 : const Element<dim>& new_element,
257 : const ::dg::MortarMap<dim, Mesh<dim>>& neighbor_mesh,
258 : const TimeStepId& current_temporal_id,
259 : const std::pair<Mesh<dim>, Element<dim>>& old_mesh_and_element) {
260 : const auto& [old_mesh, old_element] = old_mesh_and_element;
261 : ASSERT(old_element.id() == new_element.id(),
262 : "p-refinement should not have changed the element id");
263 :
264 : const bool mesh_changed = old_mesh != new_mesh;
265 :
266 : auto new_mortar_infos =
267 : detail::mortar_infos(domain, new_element, new_mesh, neighbor_mesh,
268 : Metavariables::local_time_stepping);
269 :
270 : // The old mortars must be removed from the MortarMaps before the
271 : // new ones can be added to avoid potentially exceeding the
272 : // MortarMap capacity, so we collect the new data and add it at
273 : // the end.
274 : using NewMortarEntry = std::tuple<
275 : DirectionalId<dim>, Mesh<dim - 1>,
276 : std::optional<typename mortar_data_history_type::mapped_type>>;
277 : std::vector<NewMortarEntry> new_mortars{};
278 :
279 : for (const auto& [direction, neighbors] : new_element.neighbors()) {
280 : const auto sliced_away_dimension = direction.dimension();
281 : const auto old_face_mesh = old_mesh.slice_away(sliced_away_dimension);
282 : const auto new_face_mesh = new_mesh.slice_away(sliced_away_dimension);
283 : const bool face_mesh_changed = old_face_mesh != new_face_mesh;
284 : if (face_mesh_changed) {
285 : (*normal_covector_and_magnitude)[direction] = std::nullopt;
286 : }
287 : for (const auto& neighbor : neighbors) {
288 : const DirectionalId<dim> mortar_id{direction, neighbor};
289 : const auto& new_neighbor_mesh = neighbor_mesh.at(mortar_id);
290 : const auto new_mortar_mesh = ::dg::mortar_mesh(
291 : new_face_mesh, new_neighbor_mesh.slice_away(sliced_away_dimension));
292 : if (mortar_mesh->contains(mortar_id)) {
293 : // Set the mortar mesh, but do not project any existing mesh
294 : // data. The mesh needs to have a valid value in order to
295 : // project data when we send to our neighbors. If the mesh
296 : // resolution needs to be changed afterwards because of a
297 : // change in the neighbor, that's fine, because
298 : // element-to-mortar projections are lossless. Projecting
299 : // existing data would be bad, because we might erroneously
300 : // decrease the mesh resolution, losing the high-order
301 : // modes.
302 : mortar_mesh->at(mortar_id) = new_mortar_mesh;
303 : // mortar_data does not need projecting as it has already been used
304 : // and will be resized automatically
305 : if (mesh_changed and not mortar_data_history->empty()) {
306 : auto& boundary_history = mortar_data_history->at(mortar_id);
307 : auto local_history = boundary_history.local();
308 : const auto project_local_boundary_data =
309 : [&new_face_mesh, &new_mesh](
310 : const TimeStepId& /* id */,
311 : const gsl::not_null<::evolution::dg::MortarData<dim>*>
312 : data) {
313 : return p_project_geometric_data(data, new_face_mesh,
314 : new_mesh);
315 : };
316 : local_history.for_each(project_local_boundary_data);
317 : }
318 : } else {
319 : const auto& new_mortar_size =
320 : new_mortar_infos.at(mortar_id).mortar_size();
321 :
322 : std::optional<typename mortar_data_history_type::mapped_type>
323 : new_history{};
324 : for (const auto& [old_mortar_id, old_history] :
325 : *mortar_data_history) {
326 : if (old_mortar_id.direction() != direction or
327 : not overlapping(old_mortar_id.id(), neighbor)) {
328 : continue;
329 : }
330 : const auto& old_mortar_size =
331 : mortar_infos->at(old_mortar_id).mortar_size();
332 : if (not new_history.has_value()) {
333 : new_history.emplace(old_history);
334 : new_history->remote().clear();
335 : auto local_history = new_history->local();
336 : if (mesh_changed) {
337 : const auto project_face_data =
338 : [&new_face_mesh, &new_mesh](
339 : const TimeStepId& /* id */,
340 : const gsl::not_null<::evolution::dg::MortarData<dim>*>
341 : data) {
342 : return p_project_geometric_data(data, new_face_mesh,
343 : new_mesh);
344 : };
345 : local_history.for_each(project_face_data);
346 : }
347 : const auto project_mortar_data =
348 : [&new_mortar_mesh, &new_mortar_size, &old_mortar_size](
349 : const TimeStepId& /* id */,
350 : const gsl::not_null<::evolution::dg::MortarData<dim>*>
351 : data) {
352 : const auto& old_mortar_mesh = data->mortar_mesh.value();
353 : const auto mortar_projection_matrices =
354 : Spectral::projection_matrices(
355 : old_mortar_mesh, new_mortar_mesh, old_mortar_size,
356 : new_mortar_size);
357 : DataVector& vars = data->mortar_data.value();
358 : vars = apply_matrices(mortar_projection_matrices, vars,
359 : old_mortar_mesh.extents());
360 : data->mortar_mesh = new_mortar_mesh;
361 : return true;
362 : };
363 : local_history.for_each(project_mortar_data);
364 : } else {
365 : auto local_history = new_history->local();
366 : const auto old_local_history = old_history.local();
367 : const auto project_local_mortar_data =
368 : [&new_mortar_mesh, &new_mortar_size, &old_local_history,
369 : &old_mortar_size](
370 : const TimeStepId& id,
371 : const gsl::not_null<::evolution::dg::MortarData<dim>*>
372 : data) {
373 : const auto& old_data = old_local_history.data(id);
374 : const auto& old_mortar_mesh = old_data.mortar_mesh.value();
375 : const auto mortar_projection_matrices =
376 : Spectral::projection_matrices(
377 : old_mortar_mesh, new_mortar_mesh, old_mortar_size,
378 : new_mortar_size);
379 : data->mortar_data.value() +=
380 : apply_matrices(mortar_projection_matrices,
381 : old_data.mortar_data.value(),
382 : old_mortar_mesh.extents());
383 : return true;
384 : };
385 : local_history.for_each(project_local_mortar_data);
386 : }
387 : }
388 :
389 : new_mortars.emplace_back(mortar_id, new_mortar_mesh, new_history);
390 : }
391 : }
392 : }
393 :
394 : if (not new_mortars.empty()) {
395 : domain::remove_nonexistent_neighbors(mortar_mesh, new_element);
396 : domain::remove_nonexistent_neighbors(mortar_next_temporal_id,
397 : new_element);
398 : domain::remove_nonexistent_neighbors(mortar_data_history, new_element);
399 : for (auto& [mortar_id, new_mortar_mesh, new_history] : new_mortars) {
400 : mortar_mesh->emplace(mortar_id, std::move(new_mortar_mesh));
401 : // We only do h refinement at a slab boundary, so we know all
402 : // the neighbors are aligned with us temporally.
403 : mortar_next_temporal_id->emplace(mortar_id, current_temporal_id);
404 : if (new_history.has_value()) {
405 : mortar_data_history->emplace(mortar_id, *new_history);
406 : }
407 : }
408 : }
409 :
410 : *mortar_data = detail::empty_mortar_data(new_element);
411 : *mortar_infos = std::move(new_mortar_infos);
412 :
413 : for (const auto& direction : new_element.external_boundaries()) {
414 : const auto sliced_away_dimension = direction.dimension();
415 : const auto old_face_mesh = old_mesh.slice_away(sliced_away_dimension);
416 : const auto new_face_mesh = new_mesh.slice_away(sliced_away_dimension);
417 : const bool face_mesh_changed = old_face_mesh != new_face_mesh;
418 : if (face_mesh_changed) {
419 : (*normal_covector_and_magnitude)[direction] = std::nullopt;
420 : }
421 : }
422 : }
423 :
424 : template <typename... ParentTags>
425 0 : static void apply(
426 : const gsl::not_null<
427 : ::dg::MortarMap<dim, evolution::dg::MortarDataHolder<dim>>*>
428 : mortar_data,
429 : const gsl::not_null<::dg::MortarMap<dim, Mesh<dim - 1>>*> mortar_mesh,
430 : const gsl::not_null<::dg::MortarMap<dim, MortarInfo<dim>>*> mortar_infos,
431 : const gsl::not_null<::dg::MortarMap<dim, TimeStepId>*>
432 : mortar_next_temporal_id,
433 : const gsl::not_null<
434 : DirectionMap<dim, std::optional<magnitude_and_normal_type>>*>
435 : normal_covector_and_magnitude,
436 : const gsl::not_null<mortar_data_history_type*> mortar_data_history,
437 : const Domain<dim>& domain, const Mesh<dim>& new_mesh,
438 : const Element<dim>& new_element,
439 : const ::dg::MortarMap<dim, Mesh<dim>>& neighbor_mesh,
440 : const TimeStepId& /*possibly_unset*/,
441 : const tuples::TaggedTuple<ParentTags...>& parent_items) {
442 : detail::h_refine_structure(
443 : mortar_data, mortar_mesh, mortar_infos, mortar_next_temporal_id,
444 : normal_covector_and_magnitude, domain, new_mesh, new_element,
445 : neighbor_mesh, get<::Tags::TimeStepId>(parent_items),
446 : Metavariables::local_time_stepping);
447 :
448 : const auto& old_element = get<domain::Tags::Element<dim>>(parent_items);
449 : const auto& old_histories = get<mortar_data_history_tag>(parent_items);
450 : for (const auto& [direction, neighbors] : new_element.neighbors()) {
451 : for (const auto& neighbor : neighbors) {
452 : const DirectionalId<dim> mortar_id{direction, neighbor};
453 : if (mortar_infos->at(mortar_id).time_stepping_policy() !=
454 : TimeSteppingPolicy::Conservative) {
455 : continue;
456 : }
457 : if (const auto old_history = old_histories.find(mortar_id);
458 : old_history != old_histories.end()) {
459 : // The neighbor did not h-refine, so we have to project
460 : // its mortar data from our parent.
461 : auto& new_history =
462 : mortar_data_history->emplace(mortar_id, old_history->second)
463 : .first->second;
464 : new_history.local().clear();
465 : auto remote_history = new_history.remote();
466 : const auto& new_mortar_mesh = mortar_mesh->at(mortar_id);
467 : const auto& orientation = neighbors.orientation(neighbor);
468 : const auto new_mortar_size = domain::child_size(
469 : ::dg::mortar_segments(new_element.id(), neighbor,
470 : direction.dimension(), orientation),
471 : ::dg::mortar_segments(old_element.id(), neighbor,
472 : direction.dimension(), orientation));
473 : const auto project_mortar_data =
474 : [&new_mortar_mesh, &new_mortar_size](
475 : const TimeStepId& /* id */,
476 : const gsl::not_null<::evolution::dg::MortarData<dim>*> data) {
477 : const auto& old_mortar_mesh = data->mortar_mesh.value();
478 : const auto mortar_projection_matrices =
479 : Spectral::projection_matrices(
480 : old_mortar_mesh, new_mortar_mesh,
481 : make_array<dim - 1>(Spectral::SegmentSize::Full),
482 : new_mortar_size);
483 : DataVector& vars = data->mortar_data.value();
484 : vars = apply_matrices(mortar_projection_matrices, vars,
485 : old_mortar_mesh.extents());
486 : data->mortar_mesh = new_mortar_mesh;
487 : return true;
488 : };
489 : remote_history.for_each(project_mortar_data);
490 : } else {
491 : // Neither this element nor the neighbor existed before
492 : // refinement.
493 : mortar_data_history->emplace(
494 : mortar_id, typename mortar_data_history_type::mapped_type{});
495 : }
496 : }
497 : }
498 : }
499 :
500 : template <typename... ChildTags>
501 0 : static void apply(
502 : const gsl::not_null<
503 : ::dg::MortarMap<dim, evolution::dg::MortarDataHolder<dim>>*>
504 : mortar_data,
505 : const gsl::not_null<::dg::MortarMap<dim, Mesh<dim - 1>>*> mortar_mesh,
506 : const gsl::not_null<::dg::MortarMap<dim, MortarInfo<dim>>*> mortar_infos,
507 : const gsl::not_null<::dg::MortarMap<dim, TimeStepId>*>
508 : mortar_next_temporal_id,
509 : const gsl::not_null<
510 : DirectionMap<dim, std::optional<magnitude_and_normal_type>>*>
511 : normal_covector_and_magnitude,
512 : const gsl::not_null<mortar_data_history_type*> mortar_data_history,
513 : const Domain<dim>& domain, const Mesh<dim>& new_mesh,
514 : const Element<dim>& new_element,
515 : const ::dg::MortarMap<dim, Mesh<dim>>& neighbor_mesh,
516 : const TimeStepId& /*possibly_unset*/,
517 : const std::unordered_map<
518 : ElementId<dim>, tuples::TaggedTuple<ChildTags...>>& children_items) {
519 : detail::h_refine_structure(
520 : mortar_data, mortar_mesh, mortar_infos, mortar_next_temporal_id,
521 : normal_covector_and_magnitude, domain, new_mesh, new_element,
522 : neighbor_mesh, get<::Tags::TimeStepId>(children_items.begin()->second),
523 : Metavariables::local_time_stepping);
524 :
525 : for (const auto& [direction, neighbors] : new_element.neighbors()) {
526 : for (const auto& neighbor : neighbors) {
527 : const DirectionalId<dim> mortar_id{direction, neighbor};
528 : if (mortar_infos->at(mortar_id).time_stepping_policy() !=
529 : TimeSteppingPolicy::Conservative) {
530 : continue;
531 : }
532 : std::optional<typename mortar_data_history_type::mapped_type>
533 : new_history{};
534 : for (const auto& [child, child_items] : children_items) {
535 : const auto& old_histories =
536 : get<mortar_data_history_tag>(child_items);
537 : if (const auto old_history = old_histories.find(mortar_id);
538 : old_history != old_histories.end()) {
539 : // The neighbor did not h-refine, so we have to project
540 : // its mortar data from our children.
541 : const auto& new_mortar_mesh = mortar_mesh->at(mortar_id);
542 : const auto& orientation = neighbors.orientation(neighbor);
543 : const auto old_mortar_size = domain::child_size(
544 : ::dg::mortar_segments(child, neighbor, direction.dimension(),
545 : orientation),
546 : ::dg::mortar_segments(new_element.id(), neighbor,
547 : direction.dimension(), orientation));
548 : if (not new_history.has_value()) {
549 : new_history.emplace(old_history->second);
550 : new_history->local().clear();
551 : auto remote_history = new_history->remote();
552 : const auto project_mortar_data =
553 : [&new_mortar_mesh, &old_mortar_size](
554 : const TimeStepId& /* id */,
555 : const gsl::not_null<::evolution::dg::MortarData<dim>*>
556 : data) {
557 : const auto& old_mortar_mesh = data->mortar_mesh.value();
558 : const auto mortar_projection_matrices =
559 : Spectral::projection_matrices(
560 : old_mortar_mesh, new_mortar_mesh, old_mortar_size,
561 : make_array<dim - 1>(Spectral::SegmentSize::Full));
562 : DataVector& vars = data->mortar_data.value();
563 : vars = apply_matrices(mortar_projection_matrices, vars,
564 : old_mortar_mesh.extents());
565 : data->mortar_mesh = new_mortar_mesh;
566 : return true;
567 : };
568 : remote_history.for_each(project_mortar_data);
569 : } else {
570 : auto remote_history = new_history->remote();
571 : const auto old_remote_history = old_history->second.remote();
572 : const auto project_mortar_data =
573 : [&new_mortar_mesh, &old_mortar_size, &old_remote_history](
574 : const TimeStepId& id,
575 : const gsl::not_null<::evolution::dg::MortarData<dim>*>
576 : data) {
577 : const auto& old_data = old_remote_history.data(id);
578 : const auto& old_mortar_mesh =
579 : old_data.mortar_mesh.value();
580 : const auto mortar_projection_matrices =
581 : Spectral::projection_matrices(
582 : old_mortar_mesh, new_mortar_mesh, old_mortar_size,
583 : make_array<dim - 1>(Spectral::SegmentSize::Full));
584 : data->mortar_data.value() +=
585 : apply_matrices(mortar_projection_matrices,
586 : old_data.mortar_data.value(),
587 : old_mortar_mesh.extents());
588 : return true;
589 : };
590 : remote_history.for_each(project_mortar_data);
591 : }
592 : }
593 : }
594 :
595 : if (new_history.has_value()) {
596 : mortar_data_history->emplace(mortar_id, std::move(*new_history));
597 : } else {
598 : // Neither this element nor the neighbor existed before
599 : // refinement.
600 : mortar_data_history->emplace(
601 : mortar_id, typename mortar_data_history_type::mapped_type{});
602 : }
603 : }
604 : }
605 : }
606 : };
607 : } // namespace evolution::dg::Initialization
|
