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