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