Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <optional>
9 : #include <tuple>
10 : #include <type_traits>
11 : #include <unordered_map>
12 : #include <utility>
13 : #include <vector>
14 :
15 : #include "DataStructures/DataBox/DataBox.hpp"
16 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
17 : #include "DataStructures/DataVector.hpp"
18 : #include "DataStructures/SliceVariables.hpp"
19 : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
20 : #include "DataStructures/Tensor/Tensor.hpp"
21 : #include "DataStructures/Variables.hpp"
22 : #include "DataStructures/VariablesTag.hpp"
23 : #include "Domain/Creators/Tags/Domain.hpp"
24 : #include "Domain/Creators/Tags/InitialExtents.hpp"
25 : #include "Domain/Creators/Tags/InitialRefinementLevels.hpp"
26 : #include "Domain/ElementMap.hpp"
27 : #include "Domain/FaceNormal.hpp"
28 : #include "Domain/Structure/CreateInitialMesh.hpp"
29 : #include "Domain/Structure/Direction.hpp"
30 : #include "Domain/Structure/DirectionMap.hpp"
31 : #include "Domain/Structure/Element.hpp"
32 : #include "Domain/Structure/ElementId.hpp"
33 : #include "Domain/Structure/IndexToSliceAt.hpp"
34 : #include "Domain/Tags.hpp"
35 : #include "Domain/Tags/FaceNormal.hpp"
36 : #include "Domain/Tags/Faces.hpp"
37 : #include "Domain/Tags/NeighborMesh.hpp"
38 : #include "Elliptic/DiscontinuousGalerkin/Initialization.hpp"
39 : #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/Tags.hpp"
40 : #include "Elliptic/Systems/GetFluxesComputer.hpp"
41 : #include "Elliptic/Utilities/ApplyAt.hpp"
42 : #include "Elliptic/Utilities/GetAnalyticData.hpp"
43 : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
44 : #include "NumericalAlgorithms/DiscontinuousGalerkin/ProjectToBoundary.hpp"
45 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
46 : #include "NumericalAlgorithms/Spectral/LogicalCoordinates.hpp"
47 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
48 : #include "Parallel/AlgorithmExecution.hpp"
49 : #include "ParallelAlgorithms/LinearSolver/Schwarz/OverlapHelpers.hpp"
50 : #include "ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp"
51 : #include "Utilities/Gsl.hpp"
52 : #include "Utilities/MakeArray.hpp"
53 : #include "Utilities/TMPL.hpp"
54 :
55 : /// \cond
56 : namespace Parallel {
57 : template <typename Metavariables>
58 : struct GlobalCache;
59 : } // namespace Parallel
60 : namespace tuples {
61 : template <typename... Tags>
62 : struct TaggedTuple;
63 : } // namespace tuples
64 : /// \endcond
65 :
66 : /// Actions related to the DG subdomain operator
67 1 : namespace elliptic::dg::subdomain_operator::Actions {
68 :
69 : namespace detail {
70 : // Initialize the geometry of a neighbor into which an overlap extends
71 : template <size_t Dim>
72 : struct InitializeOverlapGeometry {
73 : using return_tags =
74 : tmpl::list<elliptic::dg::subdomain_operator::Tags::ExtrudingExtent,
75 : elliptic::dg::subdomain_operator::Tags::NeighborMortars<
76 : domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>, Dim>,
77 : elliptic::dg::subdomain_operator::Tags::NeighborMortars<
78 : domain::Tags::Mesh<Dim - 1>, Dim>,
79 : elliptic::dg::subdomain_operator::Tags::NeighborMortars<
80 : ::Tags::MortarSize<Dim - 1>, Dim>>;
81 : using argument_tags =
82 : tmpl::list<domain::Tags::Element<Dim>, domain::Tags::Mesh<Dim>,
83 : domain::Tags::NeighborMesh<Dim>>;
84 : void operator()(
85 : gsl::not_null<size_t*> extruding_extent,
86 : gsl::not_null<::dg::MortarMap<Dim, Scalar<DataVector>>*>
87 : neighbor_face_normal_magnitudes,
88 : gsl::not_null<::dg::MortarMap<Dim, Mesh<Dim - 1>>*>
89 : neighbor_mortar_meshes,
90 : gsl::not_null<::dg::MortarMap<Dim, ::dg::MortarSize<Dim - 1>>*>
91 : neighbor_mortar_sizes,
92 : const Element<Dim>& element, const Mesh<Dim>& mesh,
93 : const DirectionalIdMap<Dim, Mesh<Dim>>& neighbor_meshes,
94 : const ElementId<Dim>& element_id, const Direction<Dim>& overlap_direction,
95 : std::optional<size_t> max_overlap) const;
96 : };
97 : } // namespace detail
98 :
99 : /*!
100 : * \brief Initialize the geometry for the DG subdomain operator
101 : *
102 : * Initializes tags that define the geometry of overlap regions with neighboring
103 : * elements. The data needs to be updated if the geometry of neighboring
104 : * elements changes.
105 : *
106 : * Note that the geometry depends on the system and on the choice of background
107 : * through the normalization of face normals, which involves a metric.
108 : *
109 : * The `FromInitialDomain` template parameter controls how the
110 : * next-to-nearest-neighbor information is acquired:
111 : * - If `FromInitialDomain = true`, uses the initial parameters from the input
112 : * file. This can be used at the beginning of a simulation or if the domain
113 : * never changes (no AMR).
114 : * - If `FromInitialDomain = false`, uses the mesh, element, and neighbor meshes
115 : * of each neighbor. This typically requires communication to exchange this
116 : * data between neighbors each time the domain changes, e.g. with AMR.
117 : *
118 : * DataBox:
119 : * - Uses:
120 : * - `BackgroundTag`
121 : * - `domain::Tags::Element<Dim>`
122 : * - If `FromInitialDomain = true`:
123 : * - `domain::Tags::InitialExtents<Dim>`
124 : * - `domain::Tags::InitialRefinementLevels<Dim>`
125 : * - `elliptic::dg::Tags::Quadrature`
126 : * - If `FromInitialDomain = false`:
127 : * - `Overlaps<domain::Tags::Mesh<Dim>>`
128 : * - `Overlaps<domain::Tags::Element<Dim>>`
129 : * - `Overlaps<domain::Tags::NeighborMesh<Dim>>`
130 : * - `domain::Tags::Domain<Dim>`
131 : * - `domain::Tags::FunctionsOfTime`
132 : * - `LinearSolver::Schwarz::Tags::MaxOverlap<OptionsGroup>`
133 : * - Adds: Many tags prefixed with `LinearSolver::Schwarz::Tags::Overlaps`. See
134 : * `elliptic::dg::Actions::InitializeDomain` and
135 : * `elliptic::dg::Actions::initialize_operator` for a complete list.
136 : */
137 : template <typename System, typename BackgroundTag, typename OptionsGroup,
138 : bool FromInitialDomain = true>
139 1 : struct InitializeSubdomain {
140 : private:
141 0 : static constexpr size_t Dim = System::volume_dim;
142 0 : static constexpr bool is_curved =
143 : not std::is_same_v<typename System::inv_metric_tag, void>;
144 0 : static constexpr bool has_background_fields =
145 : not std::is_same_v<typename System::background_fields, tmpl::list<>>;
146 :
147 0 : using InitializeGeometry =
148 : tmpl::conditional_t<FromInitialDomain,
149 : elliptic::dg::InitializeGeometry<Dim>,
150 : elliptic::dg::ProjectGeometry<Dim>>;
151 0 : using InitializeOverlapGeometry = detail::InitializeOverlapGeometry<Dim>;
152 0 : using InitializeFacesAndMortars = elliptic::dg::InitializeFacesAndMortars<
153 : Dim, typename System::inv_metric_tag, BackgroundTag>;
154 :
155 : template <typename Tag>
156 0 : using overlaps_tag =
157 : LinearSolver::Schwarz::Tags::Overlaps<Tag, Dim, OptionsGroup>;
158 : // Only slice those background fields to internal boundaries that are
159 : // necessary for the DG operator, i.e. the background fields in the
160 : // System::fluxes_computer::argument_tags
161 0 : using fluxes_non_background_args =
162 : tmpl::list_difference<elliptic::get_fluxes_argument_tags<System, true>,
163 : typename System::background_fields>;
164 0 : using background_fields_internal =
165 : tmpl::list_difference<elliptic::get_fluxes_argument_tags<System, true>,
166 : fluxes_non_background_args>;
167 : // Slice all background fields to external boundaries for use in boundary
168 : // conditions
169 0 : using background_fields_external = typename System::background_fields;
170 :
171 : public:
172 0 : using simple_tags_from_options =
173 : tmpl::list<domain::Tags::InitialExtents<Dim>,
174 : domain::Tags::InitialRefinementLevels<Dim>>;
175 0 : using const_global_cache_tags =
176 : tmpl::list<LinearSolver::Schwarz::Tags::MaxOverlap<OptionsGroup>>;
177 0 : using simple_tags = db::wrap_tags_in<
178 : overlaps_tag,
179 : tmpl::append<
180 : typename InitializeGeometry::return_tags,
181 : typename InitializeFacesAndMortars::return_tags,
182 : typename InitializeOverlapGeometry::return_tags,
183 : tmpl::conditional_t<
184 : has_background_fields,
185 : tmpl::list<::Tags::Variables<typename System::background_fields>>,
186 : tmpl::list<>>,
187 : domain::make_faces_tags<Dim, typename System::background_fields>>>;
188 0 : using compute_tags = tmpl::list<>;
189 :
190 : template <typename DataBox, typename... InboxTags, typename Metavariables,
191 : typename ActionList, typename ParallelComponent>
192 0 : static Parallel::iterable_action_return_t apply(
193 : DataBox& box, const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
194 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
195 : const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
196 : const ParallelComponent* const /*meta*/) {
197 : const auto& element = db::get<domain::Tags::Element<Dim>>(box);
198 : for (const auto& [direction, neighbors] : element.neighbors()) {
199 : for (const auto& neighbor_id : neighbors) {
200 : const auto& orientation = neighbors.orientation(neighbor_id);
201 : const auto direction_from_neighbor = orientation(direction.opposite());
202 : const LinearSolver::Schwarz::OverlapId<Dim> overlap_id{direction,
203 : neighbor_id};
204 : // Initialize background-agnostic geometry on overlaps
205 : elliptic::util::mutate_apply_at<
206 : db::wrap_tags_in<overlaps_tag,
207 : typename InitializeGeometry::return_tags>,
208 : tmpl::append<
209 : db::wrap_tags_in<overlaps_tag,
210 : tmpl::list_difference<
211 : typename InitializeGeometry::argument_tags,
212 : typename InitializeGeometry::volume_tags>>,
213 : typename InitializeGeometry::volume_tags>,
214 : typename InitializeGeometry::volume_tags>(
215 : InitializeGeometry{}, make_not_null(&box), overlap_id, neighbor_id);
216 : // Initialize subdomain-specific tags on overlaps
217 : elliptic::util::mutate_apply_at<
218 : db::wrap_tags_in<overlaps_tag,
219 : typename InitializeOverlapGeometry::return_tags>,
220 : db::wrap_tags_in<overlaps_tag,
221 : typename InitializeOverlapGeometry::argument_tags>,
222 : tmpl::list<>>(
223 : InitializeOverlapGeometry{}, make_not_null(&box), overlap_id,
224 : neighbor_id, direction_from_neighbor,
225 : db::get<LinearSolver::Schwarz::Tags::MaxOverlap<OptionsGroup>>(
226 : box));
227 : // Initialize faces and mortars on overlaps
228 : elliptic::util::mutate_apply_at<
229 : db::wrap_tags_in<overlaps_tag,
230 : typename InitializeFacesAndMortars::return_tags>,
231 : tmpl::append<
232 : db::wrap_tags_in<
233 : overlaps_tag,
234 : tmpl::list_difference<
235 : typename InitializeFacesAndMortars::argument_tags,
236 : typename InitializeFacesAndMortars::volume_tags>>,
237 : typename InitializeFacesAndMortars::volume_tags>,
238 : typename InitializeFacesAndMortars::volume_tags>(
239 : InitializeFacesAndMortars{}, make_not_null(&box), overlap_id);
240 : if constexpr (has_background_fields) {
241 : // Background fields
242 : initialize_background_fields(make_not_null(&box), overlap_id);
243 : }
244 : // Faces on the other side of the overlapped element's mortars
245 : initialize_remote_faces(make_not_null(&box), overlap_id);
246 : } // neighbors in direction
247 : } // directions
248 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
249 : }
250 :
251 : private:
252 : template <typename DbTagsList>
253 0 : static void initialize_background_fields(
254 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
255 : const LinearSolver::Schwarz::OverlapId<Dim>& overlap_id) {
256 : const auto& background = db::get<BackgroundTag>(*box);
257 : DirectionMap<Dim, Variables<typename System::background_fields>>
258 : face_background_fields{};
259 : elliptic::util::mutate_apply_at<
260 : tmpl::list<overlaps_tag<
261 : ::Tags::Variables<typename System::background_fields>>>,
262 : db::wrap_tags_in<
263 : overlaps_tag,
264 : tmpl::list<domain::Tags::Coordinates<Dim, Frame::Inertial>,
265 : domain::Tags::Mesh<Dim>,
266 : domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
267 : Frame::Inertial>>>,
268 : tmpl::list<>>(
269 : [&background, &face_background_fields, &box](
270 : const gsl::not_null<Variables<typename System::background_fields>*>
271 : background_fields,
272 : const tnsr::I<DataVector, Dim>& inertial_coords,
273 : const Mesh<Dim>& mesh,
274 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
275 : Frame::Inertial>& inv_jacobian) {
276 : *background_fields = elliptic::util::get_analytic_data<
277 : typename System::background_fields>(
278 : background, *box, inertial_coords, mesh, inv_jacobian);
279 : for (const auto& direction : Direction<Dim>::all_directions()) {
280 : // Slice the background fields to the face instead of evaluating
281 : // them on the face coords to avoid re-computing them, and because
282 : // this is also what the DG operator currently does. The result is
283 : // the same on Gauss-Lobatto grids, but may need adjusting when
284 : // adding support for Gauss grids.
285 : face_background_fields[direction].initialize(
286 : mesh.slice_away(direction.dimension()).number_of_grid_points());
287 : ::dg::project_contiguous_data_to_boundary(
288 : make_not_null(&face_background_fields[direction]),
289 : *background_fields, mesh, direction);
290 : }
291 : },
292 : box, overlap_id);
293 : // Move face background fields into DataBox
294 : const auto mutate_assign_face_background_field =
295 : [&box, &overlap_id, &face_background_fields](
296 : auto tag_v, const Direction<Dim>& direction) {
297 : using tag = tmpl::type_from<std::decay_t<decltype(tag_v)>>;
298 : db::mutate<overlaps_tag<domain::Tags::Faces<Dim, tag>>>(
299 : [&face_background_fields, &overlap_id,
300 : &direction](const auto stored_value) {
301 : (*stored_value)[overlap_id][direction] =
302 : get<tag>(face_background_fields.at(direction));
303 : },
304 : box);
305 : };
306 : const auto& element =
307 : db::get<overlaps_tag<domain::Tags::Element<Dim>>>(*box).at(overlap_id);
308 : for (const auto& direction : element.internal_boundaries()) {
309 : tmpl::for_each<background_fields_internal>(
310 : [&mutate_assign_face_background_field, &direction](auto tag_v) {
311 : mutate_assign_face_background_field(tag_v, direction);
312 : });
313 : }
314 : for (const auto& direction : element.external_boundaries()) {
315 : tmpl::for_each<background_fields_external>(
316 : [&mutate_assign_face_background_field, &direction](auto tag_v) {
317 : mutate_assign_face_background_field(tag_v, direction);
318 : });
319 : }
320 : }
321 :
322 : // Faces on the other side of the overlapped element's mortars
323 : template <typename DbTagsList>
324 0 : static void initialize_remote_faces(
325 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
326 : const LinearSolver::Schwarz::OverlapId<Dim>& overlap_id) {
327 : const auto& element =
328 : db::get<overlaps_tag<domain::Tags::Element<Dim>>>(*box).at(overlap_id);
329 : const auto& domain = db::get<domain::Tags::Domain<Dim>>(*box);
330 : const auto& functions_of_time =
331 : db::get<domain::Tags::FunctionsOfTime>(*box);
332 : const auto& neighbor_meshes =
333 : db::get<overlaps_tag<domain::Tags::NeighborMesh<Dim>>>(*box).at(
334 : overlap_id);
335 : for (const auto& [direction, neighbors] : element.neighbors()) {
336 : ASSERT(element.face_types().at(direction) !=
337 : domain::FaceType::MultipleNonconforming,
338 : "This code needs updating to handle nonconforming blocks");
339 : for (const auto& neighbor_id : neighbors) {
340 : const auto& orientation = neighbors.orientation(neighbor_id);
341 : const auto direction_from_neighbor = orientation(direction.opposite());
342 : const ::dg::MortarId<Dim> mortar_id{direction, neighbor_id};
343 : const auto neighbor_face_mesh =
344 : orientation(neighbor_meshes.at(mortar_id))
345 : .slice_away(direction_from_neighbor.dimension());
346 : const auto neighbor_face_logical_coords = interface_logical_coordinates(
347 : neighbor_face_mesh, direction_from_neighbor);
348 : const auto& neighbor_block = domain.blocks()[neighbor_id.block_id()];
349 : const ElementMap<Dim, Frame::Inertial> neighbor_element_map{
350 : neighbor_id, neighbor_block};
351 : const auto neighbor_face_normal = unnormalized_face_normal(
352 : neighbor_face_mesh,
353 : neighbor_element_map.inv_jacobian(neighbor_face_logical_coords, 0.,
354 : functions_of_time),
355 : direction_from_neighbor);
356 : using neighbor_face_normal_magnitudes_tag = overlaps_tag<
357 : elliptic::dg::subdomain_operator::Tags::NeighborMortars<
358 : domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>, Dim>>;
359 : if constexpr (is_curved) {
360 : const auto& background = db::get<BackgroundTag>(*box);
361 : const auto neighbor_face_inertial_coords = neighbor_element_map(
362 : neighbor_face_logical_coords, 0., functions_of_time);
363 : const auto inv_metric_on_face = get<typename System::inv_metric_tag>(
364 : elliptic::util::get_analytic_data<
365 : tmpl::list<typename System::inv_metric_tag>>(
366 : background, *box, neighbor_face_inertial_coords));
367 : elliptic::util::mutate_apply_at<
368 : tmpl::list<neighbor_face_normal_magnitudes_tag>, tmpl::list<>,
369 : tmpl::list<>>(
370 : [&neighbor_face_normal,
371 : &inv_metric_on_face](const auto neighbor_face_normal_magnitude) {
372 : magnitude(neighbor_face_normal_magnitude, neighbor_face_normal,
373 : inv_metric_on_face);
374 : },
375 : box, std::make_tuple(overlap_id, mortar_id));
376 : } else {
377 : elliptic::util::mutate_apply_at<
378 : tmpl::list<neighbor_face_normal_magnitudes_tag>, tmpl::list<>,
379 : tmpl::list<>>(
380 : [&neighbor_face_normal](
381 : const auto neighbor_face_normal_magnitude) {
382 : magnitude(neighbor_face_normal_magnitude, neighbor_face_normal);
383 : },
384 : box, std::make_tuple(overlap_id, mortar_id));
385 : }
386 : } // neighbors
387 : } // internal directions
388 : }
389 : };
390 :
391 : } // namespace elliptic::dg::subdomain_operator::Actions
|
