Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <boost/functional/hash.hpp>
7 : #include <cstddef>
8 : #include <optional>
9 : #include <ostream>
10 : #include <pup.h>
11 : #include <tuple>
12 : #include <unordered_map>
13 : #include <utility>
14 :
15 : #include "DataStructures/DataBox/DataBox.hpp"
16 : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
17 : #include "DataStructures/Variables.hpp"
18 : #include "Domain/Creators/Tags/ExternalBoundaryConditions.hpp"
19 : #include "Domain/Domain.hpp"
20 : #include "Domain/FaceNormal.hpp"
21 : #include "Domain/Structure/Direction.hpp"
22 : #include "Domain/Structure/DirectionMap.hpp"
23 : #include "Domain/Structure/ElementId.hpp"
24 : #include "Domain/Tags.hpp"
25 : #include "Domain/Tags/FaceNormal.hpp"
26 : #include "Domain/Tags/Faces.hpp"
27 : #include "Domain/Tags/SurfaceJacobian.hpp"
28 : #include "Elliptic/BoundaryConditions/ApplyBoundaryCondition.hpp"
29 : #include "Elliptic/DiscontinuousGalerkin/DgOperator.hpp"
30 : #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/Tags.hpp"
31 : #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
32 : #include "Elliptic/Systems/GetFluxesComputer.hpp"
33 : #include "Elliptic/Systems/GetSourcesComputer.hpp"
34 : #include "Elliptic/Utilities/ApplyAt.hpp"
35 : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
36 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
37 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
38 : #include "ParallelAlgorithms/LinearSolver/Schwarz/ElementCenteredSubdomainData.hpp"
39 : #include "ParallelAlgorithms/LinearSolver/Schwarz/OverlapHelpers.hpp"
40 : #include "ParallelAlgorithms/LinearSolver/Schwarz/SubdomainOperator.hpp"
41 : #include "ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp"
42 : #include "Utilities/EqualWithinRoundoff.hpp"
43 : #include "Utilities/ErrorHandling/Assert.hpp"
44 : #include "Utilities/Gsl.hpp"
45 : #include "Utilities/Serialization/PupStlCpp17.hpp"
46 : #include "Utilities/TMPL.hpp"
47 :
48 : /// Items related to the restriction of the DG operator to an element-centered
49 : /// subdomain
50 : namespace elliptic::dg::subdomain_operator {
51 :
52 : namespace detail {
53 : // Wrap the `Tag` in `LinearSolver::Schwarz::Tags::Overlaps`, except if it is
54 : // included in `TakeFromCenterTags`.
55 : template <typename Tag, typename Dim, typename OptionsGroup,
56 : typename TakeFromCenterTags>
57 : struct make_overlap_tag_impl {
58 : using type = tmpl::conditional_t<
59 : tmpl::list_contains_v<TakeFromCenterTags, Tag>, Tag,
60 : LinearSolver::Schwarz::Tags::Overlaps<Tag, Dim::value, OptionsGroup>>;
61 : };
62 :
63 : // Wrap the `Tag` in `Tags::NeighborMortars`
64 : template <typename Tag, typename Dim>
65 : struct make_neighbor_mortars_tag_impl {
66 : using type = Tags::NeighborMortars<Tag, Dim::value>;
67 : };
68 : } // namespace detail
69 :
70 : /*!
71 : * \brief The elliptic DG operator on an element-centered subdomain
72 : *
73 : * This operator is a restriction of the full (linearized) DG-operator to an
74 : * element-centered subdomain with a few points overlap into neighboring
75 : * elements. It is a `LinearSolver::Schwarz::SubdomainOperator` to be used with
76 : * the Schwarz linear solver when it solves the elliptic DG operator.
77 : *
78 : * This operator requires the following tags are available on overlap regions
79 : * with neighboring elements:
80 : *
81 : * - Geometric quantities provided by
82 : * `elliptic::dg::subdomain_operator::InitializeSubdomain`.
83 : * - All `System::fluxes_computer::argument_tags` and
84 : * `System::sources_computer::argument_tags` (or
85 : * `System::fluxes_computer_linearized::argument_tags` and
86 : * `System::sources_computer_linearized::argument_tags` for nonlinear
87 : * systems), except those listed in `ArgsTagsFromCenter`. The latter will be
88 : * taken from the central element's DataBox, so they don't need to be made
89 : * available on overlaps.
90 : * - The `System::fluxes_computer::argument_tags` (or
91 : * `System::fluxes_computer_linearized::argument_tags`) on internal and
92 : * external interfaces, except those listed in
93 : * `System::fluxes_computer::volume_tags` (or
94 : * `System::fluxes_computer_linearized::volume_tags`).
95 : *
96 : * Some of these tags may require communication between elements. For example,
97 : * nonlinear system fields are constant background fields for the linearized DG
98 : * operator, but are updated in every nonlinear solver iteration. Therefore, the
99 : * updated nonlinear fields must be communicated across overlaps between
100 : * nonlinear solver iterations. To perform the communication you can use
101 : * `LinearSolver::Schwarz::Actions::SendOverlapFields` and
102 : * `LinearSolver::Schwarz::Actions::ReceiveOverlapFields`, setting
103 : * `RestrictToOverlap` to `false`. See
104 : * `LinearSolver::Schwarz::SubdomainOperator` for details.
105 : *
106 : * \par Overriding boundary conditions
107 : * Sometimes the subdomain operator should not use the boundary conditions that
108 : * have been selected when setting up the domain. For example, when the
109 : * subdomain operator is cached between non-linear solver iterations but the
110 : * boundary conditions depend on the non-linear fields, the preconditioning can
111 : * become ineffective (see
112 : * `LinearSolver::Schwarz::Actions::ResetSubdomainSolver`). Another example is
113 : * `elliptic::subdomain_preconditioners::MinusLaplacian`, where an auxiliary
114 : * Poisson system is used for preconditioning that doesn't have boundary
115 : * conditions set up in the domain. In these cases, the boundary conditions used
116 : * for the subdomain operator can be overridden with the optional
117 : * `override_boundary_conditions` argument to the `operator()`. If the
118 : * overriding boundary conditions are different from those listed in
119 : * `Metavariables::factory_creation`, you can supply the list of
120 : * boundary-condition classes to the `BoundaryConditionClasses` template
121 : * parameter. Note that the subdomain operator always applies the _linearized_
122 : * boundary conditions.
123 : *
124 : * \warning The subdomain operator hasn't been tested with periodic boundary
125 : * conditions so far.
126 : */
127 : template <typename System, typename OptionsGroup,
128 : typename BoundaryConditionClasses = tmpl::list<>>
129 1 : struct SubdomainOperator
130 : : LinearSolver::Schwarz::SubdomainOperator<System::volume_dim> {
131 : public:
132 0 : using system = System;
133 0 : using options_group = OptionsGroup;
134 :
135 : private:
136 0 : static constexpr size_t Dim = System::volume_dim;
137 :
138 : // Operator applications happen sequentially so we don't have to keep track of
139 : // the temporal id
140 0 : static constexpr size_t temporal_id = 0;
141 :
142 : // The subdomain operator always applies the linearized DG operator
143 0 : static constexpr bool linearized = true;
144 :
145 0 : using BoundaryConditionsBase = typename System::boundary_conditions_base;
146 :
147 : // These are the arguments that we need to retrieve from the DataBox and pass
148 : // to the functions in `elliptic::dg`, both on the central element and on
149 : // neighbors
150 0 : using prepare_args_tags =
151 : tmpl::list<domain::Tags::Element<Dim>, domain::Tags::Mesh<Dim>,
152 : domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
153 : Frame::Inertial>,
154 : domain::Tags::Faces<Dim, domain::Tags::FaceNormal<Dim>>,
155 : ::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
156 : ::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>>;
157 0 : using apply_args_tags = tmpl::list<
158 : domain::Tags::Element<Dim>, domain::Tags::Mesh<Dim>,
159 : domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
160 : Frame::Inertial>,
161 : domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>,
162 : domain::Tags::DetJacobian<Frame::ElementLogical, Frame::Inertial>,
163 : domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
164 : Frame::Inertial>,
165 : domain::Tags::Faces<Dim, domain::Tags::FaceNormal<Dim>>,
166 : domain::Tags::Faces<Dim, domain::Tags::FaceNormalVector<Dim>>,
167 : domain::Tags::Faces<Dim,
168 : domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>>,
169 : domain::Tags::Faces<Dim, domain::Tags::DetSurfaceJacobian<
170 : Frame::ElementLogical, Frame::Inertial>>,
171 : domain::Tags::Faces<
172 : Dim, domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
173 : Frame::Inertial>>,
174 : ::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
175 : ::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>,
176 : ::Tags::Mortars<domain::Tags::DetSurfaceJacobian<Frame::ElementLogical,
177 : Frame::Inertial>,
178 : Dim>,
179 : ::Tags::Mortars<elliptic::dg::Tags::PenaltyFactor, Dim>,
180 : elliptic::dg::Tags::Massive, elliptic::dg::Tags::Formulation>;
181 0 : using fluxes_args_tags =
182 : elliptic::get_fluxes_argument_tags<System, linearized>;
183 0 : using sources_args_tags =
184 : elliptic::get_sources_argument_tags<System, linearized>;
185 :
186 : // We need the fluxes args also on interfaces (internal and external). The
187 : // volume tags are the subset that don't have to be taken from interfaces.
188 0 : using fluxes_args_volume_tags =
189 : elliptic::get_fluxes_volume_tags<System, linearized>;
190 :
191 : // These tags can be taken directly from the central element's DataBox, even
192 : // when evaluating neighbors
193 0 : using args_tags_from_center = tmpl::remove_duplicates<tmpl::push_back<
194 : elliptic::get_fluxes_const_global_cache_tags<System, linearized>,
195 : elliptic::get_sources_const_global_cache_tags<System, linearized>,
196 : elliptic::dg::Tags::Massive, elliptic::dg::Tags::Formulation>>;
197 :
198 : // Data on neighbors is stored in the central element's DataBox in
199 : // `LinearSolver::Schwarz::Tags::Overlaps` maps, so we wrap the argument tags
200 : // with this prefix
201 0 : using make_overlap_tag =
202 : detail::make_overlap_tag_impl<tmpl::_1, tmpl::pin<tmpl::size_t<Dim>>,
203 : tmpl::pin<OptionsGroup>,
204 : tmpl::pin<args_tags_from_center>>;
205 0 : using prepare_args_tags_overlap =
206 : tmpl::transform<prepare_args_tags, make_overlap_tag>;
207 0 : using apply_args_tags_overlap =
208 : tmpl::transform<apply_args_tags, make_overlap_tag>;
209 0 : using fluxes_args_tags_overlap =
210 : tmpl::transform<fluxes_args_tags, make_overlap_tag>;
211 0 : using sources_args_tags_overlap =
212 : tmpl::transform<sources_args_tags, make_overlap_tag>;
213 0 : using fluxes_args_tags_overlap_faces = tmpl::transform<
214 : domain::make_faces_tags<Dim, fluxes_args_tags, fluxes_args_volume_tags>,
215 : make_overlap_tag>;
216 :
217 : // We also need some data on the remote side of all neighbors' mortars. Such
218 : // data is stored in the central element's DataBox in `Tags::NeighborMortars`
219 : // maps
220 0 : using make_neighbor_mortars_tag =
221 : detail::make_neighbor_mortars_tag_impl<tmpl::_1,
222 : tmpl::pin<tmpl::size_t<Dim>>>;
223 :
224 : public:
225 : /// \warning This function is not thread-safe because it accesses mutable
226 : /// memory buffers.
227 : template <typename ResultTags, typename OperandTags, typename DbTagsList>
228 1 : void operator()(
229 : const gsl::not_null<
230 : LinearSolver::Schwarz::ElementCenteredSubdomainData<Dim, ResultTags>*>
231 : result,
232 : const LinearSolver::Schwarz::ElementCenteredSubdomainData<
233 : Dim, OperandTags>& operand,
234 : const db::DataBox<DbTagsList>& box,
235 : const std::unordered_map<std::pair<size_t, Direction<Dim>>,
236 : const BoundaryConditionsBase&,
237 : boost::hash<std::pair<size_t, Direction<Dim>>>>&
238 : override_boundary_conditions = {}) const {
239 : result->destructive_resize(operand);
240 : central_mortar_data_.clear();
241 : neighbors_mortar_data_.clear();
242 :
243 : // Used to retrieve items out of the DataBox to forward to functions. This
244 : // replaces a long series of db::get calls.
245 : const auto get_items = [](const auto&... args) {
246 : return std::forward_as_tuple(args...);
247 : };
248 :
249 : // Retrieve data out of the DataBox
250 : using tags_to_retrieve = tmpl::flatten<tmpl::list<
251 : domain::Tags::ExternalBoundaryConditions<Dim>,
252 : domain::Tags::Element<Dim>,
253 : ::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
254 : // Data on overlaps with neighbors
255 : tmpl::transform<
256 : tmpl::flatten<tmpl::list<
257 : Tags::ExtrudingExtent, domain::Tags::Element<Dim>,
258 : domain::Tags::Mesh<Dim>,
259 : domain::Tags::Faces<
260 : Dim, domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>>,
261 : ::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>,
262 : ::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>,
263 : // Data on the remote side of the neighbor's mortars
264 : tmpl::transform<
265 : tmpl::list<
266 : domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>,
267 : domain::Tags::Mesh<Dim - 1>,
268 : ::Tags::MortarSize<Dim - 1>>,
269 : make_neighbor_mortars_tag>>>,
270 : make_overlap_tag>>>;
271 : const auto& [external_boundary_conditions, central_element,
272 : central_mortar_meshes, all_overlap_extents,
273 : all_neighbor_elements, all_neighbor_meshes,
274 : all_neighbor_face_normal_magnitudes,
275 : all_neighbor_mortar_meshes, all_neighbor_mortar_sizes,
276 : all_neighbors_neighbor_face_normal_magnitudes,
277 : all_neighbors_neighbor_mortar_meshes,
278 : all_neighbors_neighbor_mortar_sizes] =
279 : db::apply<tags_to_retrieve>(get_items, box);
280 : const auto fluxes_args = db::apply<fluxes_args_tags>(get_items, box);
281 : const auto sources_args = db::apply<sources_args_tags>(get_items, box);
282 : using FluxesArgs = std::decay_t<decltype(fluxes_args)>;
283 : DirectionMap<Dim, FluxesArgs> fluxes_args_on_faces{};
284 : for (const auto& direction : Direction<Dim>::all_directions()) {
285 : fluxes_args_on_faces.emplace(
286 : direction, elliptic::util::apply_at<
287 : domain::make_faces_tags<Dim, fluxes_args_tags,
288 : fluxes_args_volume_tags>,
289 : fluxes_args_volume_tags>(get_items, box, direction));
290 : }
291 :
292 : // Setup boundary conditions
293 : const auto apply_boundary_condition =
294 : [&box, &all_boundary_conditions = external_boundary_conditions,
295 : &override_boundary_conditions](const ElementId<Dim>& local_element_id,
296 : const Direction<Dim>& local_direction,
297 : auto is_overlap, const auto& map_keys,
298 : auto&&... fields_and_fluxes) {
299 : constexpr bool is_overlap_v =
300 : std::decay_t<decltype(is_overlap)>::value;
301 : // Get boundary conditions from domain, or use overridden boundary
302 : // conditions
303 : const auto& boundary_condition = [&all_boundary_conditions,
304 : &local_element_id, &local_direction,
305 : &override_boundary_conditions]()
306 : -> const BoundaryConditionsBase& {
307 : if (not override_boundary_conditions.empty()) {
308 : const auto found_overridden_boundary_conditions =
309 : override_boundary_conditions.find(
310 : {local_element_id.block_id(), local_direction});
311 : ASSERT(found_overridden_boundary_conditions !=
312 : override_boundary_conditions.end(),
313 : "Overriding boundary conditions in subdomain operator, "
314 : "but none is available for block "
315 : << local_element_id.block_id() << " in direction "
316 : << local_direction
317 : << ". Make sure you have considered this external "
318 : "boundary of the subdomain. If this is "
319 : "intentional, add support to "
320 : "elliptic::dg::SubdomainOperator.");
321 : return found_overridden_boundary_conditions->second;
322 : }
323 : const auto& boundary_conditions =
324 : all_boundary_conditions.at(local_element_id.block_id());
325 : ASSERT(boundary_conditions.contains(local_direction),
326 : "No boundary condition is available in block "
327 : << local_element_id.block_id() << " in direction "
328 : << local_direction
329 : << ". Make sure you are setting up boundary conditions "
330 : "when creating the domain.");
331 : ASSERT(
332 : dynamic_cast<const BoundaryConditionsBase*>(
333 : boundary_conditions.at(local_direction).get()) != nullptr,
334 : "The boundary condition in block "
335 : << local_element_id.block_id() << " in direction "
336 : << local_direction
337 : << " has an unexpected type. Make sure it derives off the "
338 : "'boundary_conditions_base' class set in the system.");
339 : return dynamic_cast<const BoundaryConditionsBase&>(
340 : *boundary_conditions.at(local_direction));
341 : }();
342 : elliptic::apply_boundary_condition<
343 : linearized,
344 : tmpl::conditional_t<is_overlap_v, make_overlap_tag, void>,
345 : BoundaryConditionClasses>(
346 : boundary_condition, box, map_keys,
347 : std::forward<decltype(fields_and_fluxes)>(fields_and_fluxes)...);
348 : };
349 :
350 : // Check if the subdomain data is sparse, i.e. if some elements have zero
351 : // data. If they are, the operator is a lot cheaper to apply due to its
352 : // linearity.
353 : std::unordered_set<ElementId<Dim>> elements_in_subdomain{
354 : central_element.id()};
355 : std::unordered_set<ElementId<Dim>> elements_with_zero_data{};
356 : if (equal_within_roundoff(operand.element_data, 0.)) {
357 : elements_with_zero_data.insert(central_element.id());
358 : }
359 : for (const auto& [overlap_id, overlap_data] : operand.overlap_data) {
360 : elements_in_subdomain.insert(overlap_id.id());
361 : if (equal_within_roundoff(overlap_data, 0.)) {
362 : elements_with_zero_data.insert(overlap_id.id());
363 : }
364 : }
365 : const auto is_in_subdomain =
366 : [&elements_in_subdomain](const ElementId<Dim>& element_id) {
367 : return elements_in_subdomain.find(element_id) !=
368 : elements_in_subdomain.end();
369 : };
370 : const auto data_is_zero = [&elements_with_zero_data, &is_in_subdomain](
371 : const ElementId<Dim>& element_id) {
372 : return elements_with_zero_data.find(element_id) !=
373 : elements_with_zero_data.end() or
374 : // Data outside the subdomain is zero by definition
375 : not is_in_subdomain(element_id);
376 : };
377 : const bool central_data_is_zero = data_is_zero(central_element.id());
378 :
379 : // The subdomain operator essentially does two sweeps over all elements in
380 : // the subdomain: In the first sweep it prepares the mortar data and stores
381 : // them on both sides of all mortars, and in the second sweep it consumes
382 : // the mortar data to apply the operator. This implementation is relatively
383 : // simple because it can re-use the implementation for the parallel DG
384 : // operator. However, it is also possible to apply the subdomain operator in
385 : // a single sweep over all elements, incrementally building up the mortar
386 : // data and applying boundary corrections immediately to both adjacent
387 : // elements once the data is available. That approach is possibly a
388 : // performance optimization but requires re-implementing a lot of logic for
389 : // the DG operator here. It should be considered once the subdomain operator
390 : // has been identified as the performance bottleneck. An alternative to
391 : // optimizing the subdomain operator performance is to precondition the
392 : // subdomain solve with a _much_ simpler subdomain operator, such as a
393 : // finite-difference Laplacian, so fewer applications of the more expensive
394 : // DG subdomain operator are necessary.
395 :
396 : // 1. Prepare mortar data on all elements in the subdomain and store them on
397 : // mortars, reorienting if needed
398 : //
399 : // Prepare central element
400 : const auto apply_boundary_condition_center =
401 : [&apply_boundary_condition, &local_central_element = central_element](
402 : const Direction<Dim>& local_direction,
403 : auto&&... fields_and_fluxes) {
404 : apply_boundary_condition(
405 : local_central_element.id(), local_direction, std::false_type{},
406 : local_direction,
407 : std::forward<decltype(fields_and_fluxes)>(fields_and_fluxes)...);
408 : };
409 : db::apply<prepare_args_tags>(
410 : [this, &operand](const auto&... args) {
411 : elliptic::dg::prepare_mortar_data<System, linearized>(
412 : make_not_null(¢ral_deriv_vars_),
413 : make_not_null(¢ral_primal_fluxes_),
414 : make_not_null(¢ral_mortar_data_), operand.element_data,
415 : args...);
416 : },
417 : box, temporal_id, apply_boundary_condition_center, fluxes_args,
418 : data_is_zero);
419 : // Prepare neighbors
420 : for (const auto& [direction, neighbors] : central_element.neighbors()) {
421 : for (const auto& neighbor_id : neighbors) {
422 : const auto& orientation = neighbors.orientation(neighbor_id);
423 : const auto direction_from_neighbor = orientation(direction.opposite());
424 : const LinearSolver::Schwarz::OverlapId<Dim> overlap_id{direction,
425 : neighbor_id};
426 : const auto& overlap_extent = all_overlap_extents.at(overlap_id);
427 : const auto& neighbor = all_neighbor_elements.at(overlap_id);
428 : const auto& neighbor_mesh = all_neighbor_meshes.at(overlap_id);
429 : const auto& mortar_id = overlap_id;
430 : const auto& mortar_mesh = central_mortar_meshes.at(mortar_id);
431 : const bool neighbor_data_is_zero = data_is_zero(neighbor_id);
432 :
433 : // Intercept empty overlaps. In the unlikely case that overlaps have
434 : // zero extent, meaning no point of the neighbor is part of the
435 : // subdomain (which is fairly useless, except for testing), the
436 : // subdomain is identical to the central element and no communication
437 : // with neighbors is necessary. We can just handle the mortar between
438 : // central element and neighbor and continue.
439 : if (UNLIKELY(overlap_extent == 0)) {
440 : elliptic::dg::BoundaryData<typename System::primal_fields,
441 : typename System::primal_fluxes>
442 : remote_boundary_data{};
443 : remote_boundary_data.field_data.initialize(
444 : mortar_mesh.number_of_grid_points(), 0.);
445 : central_mortar_data_.at(mortar_id).remote_insert(
446 : temporal_id, std::move(remote_boundary_data));
447 : continue;
448 : }
449 :
450 : // Copy the central element's mortar data to the neighbor
451 : if (not(central_data_is_zero and neighbor_data_is_zero)) {
452 : auto oriented_mortar_data =
453 : central_mortar_data_.at(mortar_id).local_data(temporal_id);
454 : if (not orientation.is_aligned()) {
455 : oriented_mortar_data.orient_on_slice(
456 : mortar_mesh.extents(), direction.dimension(), orientation);
457 : }
458 : neighbors_mortar_data_[overlap_id][::dg::MortarId<Dim>{
459 : direction_from_neighbor,
460 : central_element.id()}]
461 : .remote_insert(temporal_id, std::move(oriented_mortar_data));
462 : }
463 :
464 : // Now we switch perspective to the neighbor. First, we extend the
465 : // overlap data to the full neighbor mesh by padding it with zeros. This
466 : // is necessary because spectral operators such as derivatives require
467 : // data on the full mesh.
468 : // Possible performance optimization: this copy can be avoided if the
469 : // overlap covers the full neighbor, so no padding with zeros is needed.
470 : if (not neighbor_data_is_zero) {
471 : LinearSolver::Schwarz::extended_overlap_data(
472 : make_not_null(&extended_operand_vars_[overlap_id]),
473 : operand.overlap_data.at(overlap_id), neighbor_mesh.extents(),
474 : overlap_extent, direction_from_neighbor);
475 : }
476 :
477 : const auto apply_boundary_condition_neighbor =
478 : [&apply_boundary_condition, &local_neighbor_id = neighbor_id,
479 : &overlap_id](const Direction<Dim>& local_direction,
480 : auto&&... fields_and_fluxes) {
481 : apply_boundary_condition(
482 : local_neighbor_id, local_direction, std::true_type{},
483 : std::forward_as_tuple(overlap_id, local_direction),
484 : std::forward<decltype(fields_and_fluxes)>(
485 : fields_and_fluxes)...);
486 : };
487 :
488 : const auto fluxes_args_on_overlap =
489 : elliptic::util::apply_at<fluxes_args_tags_overlap,
490 : args_tags_from_center>(get_items, box,
491 : overlap_id);
492 :
493 : elliptic::util::apply_at<prepare_args_tags_overlap,
494 : args_tags_from_center>(
495 : [this, &overlap_id](const auto&... args) {
496 : elliptic::dg::prepare_mortar_data<System, linearized>(
497 : make_not_null(&neighbors_deriv_vars_[overlap_id]),
498 : make_not_null(&neighbors_primal_fluxes_[overlap_id]),
499 : make_not_null(&neighbors_mortar_data_[overlap_id]),
500 : extended_operand_vars_[overlap_id], args...);
501 : },
502 : box, overlap_id, temporal_id, apply_boundary_condition_neighbor,
503 : fluxes_args_on_overlap, data_is_zero);
504 :
505 : // Copy this neighbor's mortar data to the other side of the mortars. On
506 : // the other side we either have the central element, or another element
507 : // that may or may not be part of the subdomain.
508 : const auto& neighbor_mortar_meshes =
509 : all_neighbor_mortar_meshes.at(overlap_id);
510 : for (const auto& neighbor_mortar_id_and_data :
511 : neighbors_mortar_data_.at(overlap_id)) {
512 : // No structured bindings because capturing these in lambdas doesn't
513 : // work until C++20
514 : const auto& neighbor_mortar_id = neighbor_mortar_id_and_data.first;
515 : const auto& neighbor_mortar_data = neighbor_mortar_id_and_data.second;
516 : const auto& neighbor_direction = neighbor_mortar_id.direction();
517 : const auto& neighbors_neighbor_id = neighbor_mortar_id.id();
518 : // No need to do anything on external boundaries
519 : if (neighbors_neighbor_id == ElementId<Dim>::external_boundary_id()) {
520 : continue;
521 : }
522 : const auto& neighbor_orientation =
523 : neighbor.neighbors()
524 : .at(neighbor_direction)
525 : .orientation(neighbors_neighbor_id);
526 : const auto neighbors_neighbor_direction =
527 : neighbor_orientation(neighbor_direction.opposite());
528 : const ::dg::MortarId<Dim> mortar_id_from_neighbors_neighbor{
529 : neighbors_neighbor_direction, neighbor_id};
530 : const auto send_mortar_data =
531 : [&neighbor_orientation, &neighbor_mortar_meshes,
532 : &neighbor_mortar_data, &neighbor_mortar_id, &neighbor_direction,
533 : &neighbor_data_is_zero,
534 : &data_is_zero](auto& remote_mortar_data,
535 : const ElementId<Dim>& remote_element_id) {
536 : if (neighbor_data_is_zero and data_is_zero(remote_element_id)) {
537 : return;
538 : }
539 : const auto& neighbor_mortar_mesh =
540 : neighbor_mortar_meshes.at(neighbor_mortar_id);
541 : auto oriented_neighbor_mortar_data =
542 : neighbor_mortar_data.local_data(temporal_id);
543 : if (not neighbor_orientation.is_aligned()) {
544 : oriented_neighbor_mortar_data.orient_on_slice(
545 : neighbor_mortar_mesh.extents(),
546 : neighbor_direction.dimension(), neighbor_orientation);
547 : }
548 : remote_mortar_data.remote_insert(
549 : temporal_id, std::move(oriented_neighbor_mortar_data));
550 : };
551 : if (neighbors_neighbor_id == central_element.id() and
552 : mortar_id_from_neighbors_neighbor == mortar_id) {
553 : send_mortar_data(central_mortar_data_.at(mortar_id),
554 : central_element.id());
555 : continue;
556 : }
557 : // Determine whether the neighbor's neighbor overlaps with the
558 : // subdomain and find its overlap ID if it does.
559 : const auto neighbors_neighbor_overlap_id =
560 : [&local_all_neighbor_mortar_meshes = all_neighbor_mortar_meshes,
561 : &neighbors_neighbor_id, &mortar_id_from_neighbors_neighbor,
562 : &is_in_subdomain]()
563 : -> std::optional<LinearSolver::Schwarz::OverlapId<Dim>> {
564 : if (not is_in_subdomain(neighbors_neighbor_id)) {
565 : return std::nullopt;
566 : }
567 : for (const auto& [local_overlap_id, local_mortar_meshes] :
568 : local_all_neighbor_mortar_meshes) {
569 : if (local_overlap_id.id() != neighbors_neighbor_id) {
570 : continue;
571 : }
572 : for (const auto& local_mortar_id_and_mesh : local_mortar_meshes) {
573 : if (local_mortar_id_and_mesh.first ==
574 : mortar_id_from_neighbors_neighbor) {
575 : return local_overlap_id;
576 : }
577 : }
578 : }
579 : ERROR("The neighbor's neighbor "
580 : << neighbors_neighbor_id
581 : << " is part of the subdomain, but we didn't find its "
582 : "overlap ID. This is a bug, so please file an issue.");
583 : }();
584 : if (neighbors_neighbor_overlap_id.has_value()) {
585 : // The neighbor's neighbor is part of the subdomain so we copy the
586 : // mortar data over. Once the loop is complete we will also have
587 : // received mortar data back. At that point, both neighbors have a
588 : // copy of each other's mortar data, which is the subject of the
589 : // possible optimizations mentioned above. Note that the data may
590 : // differ by orientations.
591 : send_mortar_data(
592 : neighbors_mortar_data_[*neighbors_neighbor_overlap_id]
593 : [mortar_id_from_neighbors_neighbor],
594 : neighbors_neighbor_overlap_id->id());
595 : } else if (not neighbor_data_is_zero) {
596 : // The neighbor's neighbor does not overlap with the subdomain, so
597 : // we don't copy mortar data and also don't expect to receive any.
598 : // Instead, we assume the data on it is zero and manufacture
599 : // appropriate remote boundary data.
600 : const auto& neighbors_neighbor_mortar_mesh =
601 : all_neighbors_neighbor_mortar_meshes.at(overlap_id)
602 : .at(neighbor_mortar_id);
603 : elliptic::dg::BoundaryData<typename System::primal_fields,
604 : typename System::primal_fluxes>
605 : zero_mortar_data{};
606 : zero_mortar_data.field_data.initialize(
607 : neighbors_neighbor_mortar_mesh.number_of_grid_points(), 0.);
608 : neighbors_mortar_data_.at(overlap_id)
609 : .at(neighbor_mortar_id)
610 : .remote_insert(temporal_id, std::move(zero_mortar_data));
611 : }
612 : } // loop over neighbor's mortars
613 : } // loop over neighbors
614 : } // loop over directions
615 :
616 : // 2. Apply the operator on all elements in the subdomain
617 : //
618 : // Apply on central element
619 : db::apply<apply_args_tags>(
620 : [this, &result, &operand](const auto&... args) {
621 : elliptic::dg::apply_operator<System, linearized>(
622 : make_not_null(&result->element_data),
623 : make_not_null(¢ral_mortar_data_), operand.element_data,
624 : central_primal_fluxes_, args...);
625 : },
626 : box, temporal_id, fluxes_args_on_faces, sources_args, data_is_zero);
627 : // Apply on neighbors
628 : for (const auto& [direction, neighbors] : central_element.neighbors()) {
629 : for (const auto& neighbor_id : neighbors) {
630 : const auto& orientation = neighbors.orientation(neighbor_id);
631 : const auto direction_from_neighbor = orientation(direction.opposite());
632 : const LinearSolver::Schwarz::OverlapId<Dim> overlap_id{direction,
633 : neighbor_id};
634 : const auto& overlap_extent = all_overlap_extents.at(overlap_id);
635 : const auto& neighbor_mesh = all_neighbor_meshes.at(overlap_id);
636 :
637 : if (UNLIKELY(overlap_extent == 0)) {
638 : continue;
639 : }
640 :
641 : DirectionMap<Dim, FluxesArgs> fluxes_args_on_overlap_faces{};
642 : for (const auto& neighbor_direction :
643 : Direction<Dim>::all_directions()) {
644 : fluxes_args_on_overlap_faces.emplace(
645 : neighbor_direction,
646 : elliptic::util::apply_at<fluxes_args_tags_overlap_faces,
647 : args_tags_from_center>(
648 : get_items, box,
649 : std::forward_as_tuple(overlap_id, neighbor_direction)));
650 : }
651 :
652 : elliptic::util::apply_at<apply_args_tags_overlap,
653 : args_tags_from_center>(
654 : [this, &overlap_id](const auto&... args) {
655 : elliptic::dg::apply_operator<System, linearized>(
656 : make_not_null(&extended_results_[overlap_id]),
657 : make_not_null(&neighbors_mortar_data_.at(overlap_id)),
658 : extended_operand_vars_.at(overlap_id),
659 : neighbors_primal_fluxes_.at(overlap_id), args...);
660 : },
661 : box, overlap_id, temporal_id, fluxes_args_on_overlap_faces,
662 : elliptic::util::apply_at<sources_args_tags_overlap,
663 : args_tags_from_center>(get_items, box,
664 : overlap_id),
665 : data_is_zero);
666 :
667 : // Restrict the extended operator data back to the subdomain, assuming
668 : // we can discard any data outside the overlaps. WARNING: This
669 : // assumption may break with changes to the DG operator that affect its
670 : // sparsity. For example, multiplying the DG operator with the _full_
671 : // inverse mass-matrix ("massless" scheme with no "mass-lumping"
672 : // approximation) means that lifted boundary corrections bleed into the
673 : // volume.
674 : // Possible performance optimization: this copy can be avoided if the
675 : // overlap covers the full neighbor, so no restriction is needed.
676 : if (UNLIKELY(
677 : result->overlap_data[overlap_id].number_of_grid_points() !=
678 : operand.overlap_data.at(overlap_id).number_of_grid_points())) {
679 : result->overlap_data[overlap_id].initialize(
680 : operand.overlap_data.at(overlap_id).number_of_grid_points());
681 : }
682 : LinearSolver::Schwarz::data_on_overlap(
683 : make_not_null(&result->overlap_data[overlap_id]),
684 : extended_results_.at(overlap_id), neighbor_mesh.extents(),
685 : overlap_extent, direction_from_neighbor);
686 : } // loop over neighbors
687 : } // loop over directions
688 : }
689 :
690 : // NOLINTNEXTLINE(google-runtime-references)
691 0 : void pup(PUP::er& /*p*/) {}
692 :
693 : private:
694 : // Memory buffers for repeated operator applications
695 : // NOLINTNEXTLINE(spectre-mutable)
696 : mutable Variables<
697 : db::wrap_tags_in<::Tags::deriv, typename System::primal_fields,
698 : tmpl::size_t<Dim>, Frame::Inertial>>
699 0 : central_deriv_vars_{};
700 : // NOLINTNEXTLINE(spectre-mutable)
701 0 : mutable Variables<typename System::primal_fluxes> central_primal_fluxes_{};
702 : // NOLINTNEXTLINE(spectre-mutable)
703 : mutable LinearSolver::Schwarz::OverlapMap<
704 : Dim,
705 : Variables<db::wrap_tags_in<::Tags::deriv, typename System::primal_fields,
706 : tmpl::size_t<Dim>, Frame::Inertial>>>
707 0 : neighbors_deriv_vars_{};
708 : // NOLINTNEXTLINE(spectre-mutable)
709 : mutable LinearSolver::Schwarz::OverlapMap<
710 : Dim, Variables<typename System::primal_fluxes>>
711 0 : neighbors_primal_fluxes_{};
712 : // NOLINTNEXTLINE(spectre-mutable)
713 : mutable LinearSolver::Schwarz::OverlapMap<
714 : Dim, Variables<typename System::primal_fields>>
715 0 : extended_operand_vars_{};
716 : // NOLINTNEXTLINE(spectre-mutable)
717 : mutable LinearSolver::Schwarz::OverlapMap<
718 : Dim, Variables<typename System::primal_fields>>
719 0 : extended_results_{};
720 : // NOLINTNEXTLINE(spectre-mutable)
721 : mutable ::dg::MortarMap<
722 : Dim, elliptic::dg::MortarData<size_t, typename System::primal_fields,
723 : typename System::primal_fluxes>>
724 0 : central_mortar_data_{};
725 : // NOLINTNEXTLINE(spectre-mutable)
726 : mutable LinearSolver::Schwarz::OverlapMap<
727 : Dim, ::dg::MortarMap<Dim, elliptic::dg::MortarData<
728 : size_t, typename System::primal_fields,
729 : typename System::primal_fluxes>>>
730 0 : neighbors_mortar_data_{};
731 : };
732 :
733 : } // namespace elliptic::dg::subdomain_operator
|
