Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <optional>
8 : #include <string>
9 : #include <tuple>
10 : #include <type_traits>
11 : #include <utility>
12 :
13 : #include "DataStructures/DataBox/DataBox.hpp"
14 : #include "DataStructures/DataBox/DataBoxTag.hpp"
15 : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
16 : #include "Domain/Creators/Tags/ExternalBoundaryConditions.hpp"
17 : #include "Domain/Creators/Tags/InitialExtents.hpp"
18 : #include "Domain/FaceNormal.hpp"
19 : #include "Domain/Structure/Direction.hpp"
20 : #include "Domain/Structure/DirectionMap.hpp"
21 : #include "Domain/Structure/DirectionalIdMap.hpp"
22 : #include "Domain/Structure/Element.hpp"
23 : #include "Domain/Tags.hpp"
24 : #include "Domain/Tags/FaceNormal.hpp"
25 : #include "Domain/Tags/Faces.hpp"
26 : #include "Domain/Tags/SurfaceJacobian.hpp"
27 : #include "Elliptic/BoundaryConditions/ApplyBoundaryCondition.hpp"
28 : #include "Elliptic/DiscontinuousGalerkin/DgOperator.hpp"
29 : #include "Elliptic/DiscontinuousGalerkin/Initialization.hpp"
30 : #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
31 : #include "Elliptic/Systems/GetFluxesComputer.hpp"
32 : #include "Elliptic/Systems/GetModifyBoundaryData.hpp"
33 : #include "Elliptic/Systems/GetSourcesComputer.hpp"
34 : #include "Elliptic/Utilities/ApplyAt.hpp"
35 : #include "NumericalAlgorithms/Convergence/Tags.hpp"
36 : #include "NumericalAlgorithms/DiscontinuousGalerkin/HasReceivedFromAllMortars.hpp"
37 : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
38 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
39 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
40 : #include "Parallel/AlgorithmExecution.hpp"
41 : #include "Parallel/GlobalCache.hpp"
42 : #include "Parallel/InboxInserters.hpp"
43 : #include "Parallel/Invoke.hpp"
44 : #include "ParallelAlgorithms/Amr/Projectors/CopyFromCreatorOrLeaveAsIs.hpp"
45 : #include "ParallelAlgorithms/Amr/Projectors/DefaultInitialize.hpp"
46 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
47 : #include "Utilities/ErrorHandling/Assert.hpp"
48 : #include "Utilities/GetOutput.hpp"
49 : #include "Utilities/Gsl.hpp"
50 : #include "Utilities/Literals.hpp"
51 : #include "Utilities/TMPL.hpp"
52 : #include "Utilities/TaggedTuple.hpp"
53 :
54 : /// Actions related to elliptic discontinuous Galerkin schemes
55 1 : namespace elliptic::dg::Actions {
56 : // The individual actions in this namespace are not exposed publicly because
57 : // they don't work on their own. Instead, the public interface (defined below)
58 : // exposes them in action lists.
59 : namespace detail {
60 :
61 : // This is a global ID that identifies the DG operator application. It
62 : // increments with each application of the operator.
63 : struct DgOperatorLabel {};
64 : using TemporalIdTag = Convergence::Tags::IterationId<DgOperatorLabel>;
65 :
66 : // This tag is used to communicate mortar data across neighbors
67 : template <size_t Dim, typename TemporalIdTag, typename PrimalFields,
68 : typename PrimalFluxes>
69 : struct MortarDataInboxTag
70 : : public Parallel::InboxInserters::Map<
71 : MortarDataInboxTag<Dim, TemporalIdTag, PrimalFields, PrimalFluxes>> {
72 : using temporal_id = typename TemporalIdTag::type;
73 : using type =
74 : std::map<temporal_id,
75 : DirectionalIdMap<Dim, elliptic::dg::BoundaryData<PrimalFields,
76 : PrimalFluxes>>>;
77 : };
78 :
79 : // Initializes all quantities the DG operator needs on internal and external
80 : // faces, as well as the mortars between neighboring elements. Also initializes
81 : // the variable-independent background fields in the PDEs.
82 : template <typename System, typename BackgroundTag>
83 : struct InitializeFacesMortarsAndBackground {
84 : private:
85 : static constexpr size_t Dim = System::volume_dim;
86 : static constexpr bool has_background_fields =
87 : not std::is_same_v<typename System::background_fields, tmpl::list<>>;
88 : static_assert(
89 : not(has_background_fields and std::is_same_v<BackgroundTag, void>),
90 : "The system has background fields that need initialization. Supply a "
91 : "'BackgroundTag' to 'elliptic::dg::Actions::initialize_operator'.");
92 :
93 : using InitializeFacesAndMortars = elliptic::dg::InitializeFacesAndMortars<
94 : Dim, typename System::inv_metric_tag, BackgroundTag>;
95 : using InitializeBackground = elliptic::dg::InitializeBackground<
96 : Dim, typename System::background_fields, BackgroundTag>;
97 :
98 : public:
99 : using simple_tags = tmpl::append<
100 : typename InitializeFacesAndMortars::return_tags,
101 : tmpl::conditional_t<has_background_fields,
102 : typename InitializeBackground::return_tags,
103 : tmpl::list<>>>;
104 : using compute_tags = tmpl::list<>;
105 : using const_global_cache_tags =
106 : tmpl::conditional_t<has_background_fields, tmpl::list<BackgroundTag>,
107 : tmpl::list<>>;
108 :
109 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
110 : typename ArrayIndex, typename ActionList,
111 : typename ParallelComponent>
112 : static Parallel::iterable_action_return_t apply(
113 : db::DataBox<DbTagsList>& box,
114 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
115 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
116 : const ArrayIndex& /*array_index*/, ActionList /*meta*/,
117 : const ParallelComponent* const /*meta*/) {
118 : // Initialize faces and mortars
119 : db::mutate_apply<InitializeFacesAndMortars>(make_not_null(&box));
120 : if constexpr (has_background_fields) {
121 : // Initialize background fields
122 : db::mutate_apply<InitializeBackground>(make_not_null(&box));
123 : }
124 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
125 : }
126 : };
127 :
128 : // Compute auxiliary variables and fluxes from the primal variables, prepare the
129 : // local side of all mortars and send the local mortar data to neighbors. Also
130 : // handle boundary conditions by preparing the exterior ("ghost") side of
131 : // external mortars.
132 : template <typename System, bool Linearized, typename TemporalIdTag,
133 : typename PrimalFieldsTag, typename PrimalFluxesTag,
134 : typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
135 : typename PrimalMortarFluxesTag,
136 : typename FluxesArgsTags =
137 : elliptic::get_fluxes_argument_tags<System, Linearized>,
138 : typename SourcesArgsTags =
139 : elliptic::get_sources_argument_tags<System, Linearized>>
140 : struct PrepareAndSendMortarData;
141 :
142 : template <typename System, bool Linearized, typename TemporalIdTag,
143 : typename PrimalFieldsTag, typename PrimalFluxesTag,
144 : typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
145 : typename PrimalMortarFluxesTag, typename... FluxesArgsTags,
146 : typename... SourcesArgsTags>
147 : struct PrepareAndSendMortarData<
148 : System, Linearized, TemporalIdTag, PrimalFieldsTag, PrimalFluxesTag,
149 : OperatorAppliedToFieldsTag, PrimalMortarFieldsTag, PrimalMortarFluxesTag,
150 : tmpl::list<FluxesArgsTags...>, tmpl::list<SourcesArgsTags...>> {
151 : private:
152 : static constexpr size_t Dim = System::volume_dim;
153 : using all_mortar_data_tag = ::Tags::Mortars<
154 : elliptic::dg::Tags::MortarData<typename TemporalIdTag::type,
155 : typename PrimalMortarFieldsTag::tags_list,
156 : typename PrimalMortarFluxesTag::tags_list>,
157 : Dim>;
158 : using mortar_data_inbox_tag =
159 : MortarDataInboxTag<Dim, TemporalIdTag,
160 : typename PrimalMortarFieldsTag::tags_list,
161 : typename PrimalMortarFluxesTag::tags_list>;
162 : using BoundaryConditionsBase = typename System::boundary_conditions_base;
163 : using DerivFieldsTag = db::add_tag_prefix<::Tags::deriv, PrimalFieldsTag,
164 : tmpl::size_t<Dim>, Frame::Inertial>;
165 :
166 : public:
167 : // Request these tags be added to the DataBox. We
168 : // don't actually need to initialize them, because the
169 : // `PrimalFieldsTag` will be set by other actions before applying the operator
170 : // and the remaining tags hold output of the operator.
171 : using simple_tags = tmpl::list<TemporalIdTag, PrimalFieldsTag, DerivFieldsTag,
172 : PrimalFluxesTag, OperatorAppliedToFieldsTag,
173 : all_mortar_data_tag>;
174 : using compute_tags = tmpl::list<>;
175 : using const_global_cache_tags =
176 : tmpl::list<domain::Tags::ExternalBoundaryConditions<Dim>>;
177 :
178 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
179 : typename ActionList, typename ParallelComponent>
180 : static Parallel::iterable_action_return_t apply(
181 : db::DataBox<DbTagsList>& box,
182 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
183 : Parallel::GlobalCache<Metavariables>& cache,
184 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
185 : const ParallelComponent* const /*meta*/) {
186 : const auto& temporal_id = db::get<TemporalIdTag>(box);
187 : const auto& element = db::get<domain::Tags::Element<Dim>>(box);
188 :
189 : const auto& mortar_meshes =
190 : db::get<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>>(box);
191 : const auto& boundary_conditions =
192 : db::get<domain::Tags::ExternalBoundaryConditions<Dim>>(box).at(
193 : element_id.block_id());
194 : const auto apply_boundary_condition = [&box, &boundary_conditions,
195 : &element_id](
196 : const Direction<Dim>& direction,
197 : auto&&... fields_and_fluxes) {
198 : ASSERT(boundary_conditions.contains(direction),
199 : "No boundary condition is available in block "
200 : << element_id.block_id() << " in direction " << direction
201 : << ". Make sure you are setting up boundary conditions when "
202 : "creating the domain.");
203 : ASSERT(dynamic_cast<const BoundaryConditionsBase*>(
204 : boundary_conditions.at(direction).get()) != nullptr,
205 : "The boundary condition in block "
206 : << element_id.block_id() << " in direction " << direction
207 : << " has an unexpected type. Make sure it derives off the "
208 : "'boundary_conditions_base' class set in the system.");
209 : const auto& boundary_condition =
210 : dynamic_cast<const BoundaryConditionsBase&>(
211 : *boundary_conditions.at(direction));
212 : elliptic::apply_boundary_condition<Linearized>(
213 : boundary_condition, box, direction,
214 : std::forward<decltype(fields_and_fluxes)>(fields_and_fluxes)...);
215 : };
216 :
217 : // Can't `db::get` the arguments for the boundary conditions within
218 : // `db::mutate`, so we extract the data to mutate and move it back in when
219 : // we're done.
220 : // Possible optimization: also keep memory for the mortar data around in the
221 : // DataBox. Currently the mortar data is extracted and erased by
222 : // `apply_operator` anyway, so we can just create a new map here to avoid
223 : // dealing with AMR resizing the mortars. When we keep the memory around,
224 : // its size has to be adjusted when the mesh changes during AMR.
225 : typename PrimalFluxesTag::type primal_fluxes;
226 : typename DerivFieldsTag::type deriv_fields;
227 : typename all_mortar_data_tag::type all_mortar_data{};
228 : db::mutate<PrimalFluxesTag, DerivFieldsTag>(
229 : [&primal_fluxes, &deriv_fields](const auto local_primal_fluxes,
230 : const auto local_deriv_fields) {
231 : primal_fluxes = std::move(*local_primal_fluxes);
232 : deriv_fields = std::move(*local_deriv_fields);
233 : },
234 : make_not_null(&box));
235 :
236 : elliptic::dg::prepare_mortar_data<System, Linearized>(
237 : make_not_null(&deriv_fields), make_not_null(&primal_fluxes),
238 : make_not_null(&all_mortar_data), db::get<PrimalFieldsTag>(box), element,
239 : db::get<domain::Tags::Mesh<Dim>>(box),
240 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
241 : Frame::Inertial>>(box),
242 : db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormal<Dim>>>(box),
243 : mortar_meshes,
244 : db::get<::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>>(box),
245 : temporal_id, apply_boundary_condition,
246 : std::forward_as_tuple(db::get<FluxesArgsTags>(box)...));
247 :
248 : // Move the mutated data back into the DataBox
249 : db::mutate<PrimalFluxesTag, DerivFieldsTag, all_mortar_data_tag>(
250 : [&primal_fluxes, &deriv_fields, &all_mortar_data](
251 : const auto local_primal_fluxes, const auto local_deriv_fields,
252 : const auto local_all_mortar_data) {
253 : *local_primal_fluxes = std::move(primal_fluxes);
254 : *local_deriv_fields = std::move(deriv_fields);
255 : *local_all_mortar_data = std::move(all_mortar_data);
256 : },
257 : make_not_null(&box));
258 :
259 : // Send mortar data to neighbors
260 : auto& receiver_proxy =
261 : Parallel::get_parallel_component<ParallelComponent>(cache);
262 : for (const auto& [direction, neighbors] : element.neighbors()) {
263 : const size_t dimension = direction.dimension();
264 : for (const auto& neighbor_id : neighbors) {
265 : const auto& orientation = neighbors.orientation(neighbor_id);
266 : const auto direction_from_neighbor = orientation(direction.opposite());
267 : const ::dg::MortarId<Dim> mortar_id{direction, neighbor_id};
268 : // Make a copy of the local boundary data on the mortar to send to the
269 : // neighbor
270 : auto remote_boundary_data_on_mortar =
271 : get<all_mortar_data_tag>(box).at(mortar_id).local_data(temporal_id);
272 : // Reorient the data to the neighbor orientation if necessary
273 : if (not orientation.is_aligned()) {
274 : remote_boundary_data_on_mortar.orient_on_slice(
275 : mortar_meshes.at(mortar_id).extents(), dimension, orientation);
276 : }
277 : // Send remote data to neighbor
278 : Parallel::receive_data<mortar_data_inbox_tag>(
279 : receiver_proxy[neighbor_id], temporal_id,
280 : std::make_pair(
281 : ::dg::MortarId<Dim>{direction_from_neighbor, element.id()},
282 : std::move(remote_boundary_data_on_mortar)));
283 : } // loop over neighbors in direction
284 : } // loop over directions
285 :
286 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
287 : }
288 : };
289 :
290 : // Wait until all mortar data from neighbors is available. Then add boundary
291 : // corrections to the primal fluxes, compute their derivatives (i.e. the second
292 : // derivatives of the primal variables) and add boundary corrections to the
293 : // result.
294 : template <typename System, bool Linearized, typename TemporalIdTag,
295 : typename PrimalFieldsTag, typename PrimalFluxesTag,
296 : typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
297 : typename PrimalMortarFluxesTag,
298 : typename FluxesArgsTags =
299 : elliptic::get_fluxes_argument_tags<System, Linearized>,
300 : typename SourcesArgsTags =
301 : elliptic::get_sources_argument_tags<System, Linearized>>
302 : struct ReceiveMortarDataAndApplyOperator;
303 :
304 : template <typename System, bool Linearized, typename TemporalIdTag,
305 : typename PrimalFieldsTag, typename PrimalFluxesTag,
306 : typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
307 : typename PrimalMortarFluxesTag, typename... FluxesArgsTags,
308 : typename... SourcesArgsTags>
309 : struct ReceiveMortarDataAndApplyOperator<
310 : System, Linearized, TemporalIdTag, PrimalFieldsTag, PrimalFluxesTag,
311 : OperatorAppliedToFieldsTag, PrimalMortarFieldsTag, PrimalMortarFluxesTag,
312 : tmpl::list<FluxesArgsTags...>, tmpl::list<SourcesArgsTags...>> {
313 : private:
314 : static constexpr size_t Dim = System::volume_dim;
315 : using all_mortar_data_tag = ::Tags::Mortars<
316 : elliptic::dg::Tags::MortarData<typename TemporalIdTag::type,
317 : typename PrimalMortarFieldsTag::tags_list,
318 : typename PrimalMortarFluxesTag::tags_list>,
319 : Dim>;
320 : using mortar_data_inbox_tag =
321 : MortarDataInboxTag<Dim, TemporalIdTag,
322 : typename PrimalMortarFieldsTag::tags_list,
323 : typename PrimalMortarFluxesTag::tags_list>;
324 : using DerivFieldsTag = db::add_tag_prefix<::Tags::deriv, PrimalFieldsTag,
325 : tmpl::size_t<Dim>, Frame::Inertial>;
326 : public:
327 : using const_global_cache_tags =
328 : tmpl::list<elliptic::dg::Tags::PenaltyParameter,
329 : elliptic::dg::Tags::Massive, elliptic::dg::Tags::Formulation>;
330 : using inbox_tags = tmpl::list<mortar_data_inbox_tag>;
331 :
332 : template <typename DbTags, typename... InboxTags, typename Metavariables,
333 : typename ArrayIndex, typename ActionList,
334 : typename ParallelComponent>
335 : static Parallel::iterable_action_return_t apply(
336 : db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
337 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
338 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
339 : const ParallelComponent* const /*meta*/) {
340 : const auto& temporal_id = get<TemporalIdTag>(box);
341 : const auto& element = get<domain::Tags::Element<Dim>>(box);
342 :
343 : if (not::dg::has_received_from_all_mortars<mortar_data_inbox_tag>(
344 : temporal_id, element, inboxes)) {
345 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
346 : }
347 :
348 : // Move received "remote" mortar data into the DataBox
349 : if (LIKELY(element.number_of_neighbors() > 0)) {
350 : auto received_mortar_data =
351 : std::move(tuples::get<mortar_data_inbox_tag>(inboxes)
352 : .extract(temporal_id)
353 : .mapped());
354 : db::mutate<all_mortar_data_tag>(
355 : [&received_mortar_data, &temporal_id](const auto all_mortar_data) {
356 : for (auto& [mortar_id, mortar_data] : received_mortar_data) {
357 : all_mortar_data->at(mortar_id).remote_insert(
358 : temporal_id, std::move(mortar_data));
359 : }
360 : },
361 : make_not_null(&box));
362 : }
363 :
364 : // Used to retrieve items out of the DataBox to forward to functions
365 : const auto get_items = [](const auto&... args) {
366 : return std::forward_as_tuple(args...);
367 : };
368 :
369 : // Apply DG operator
370 : using fluxes_args_tags =
371 : typename elliptic::get_fluxes_argument_tags<System, Linearized>;
372 : using fluxes_args_volume_tags =
373 : typename elliptic::get_fluxes_volume_tags<System, Linearized>;
374 : DirectionMap<Dim, std::tuple<decltype(db::get<FluxesArgsTags>(box))...>>
375 : fluxes_args_on_faces{};
376 : for (const auto& direction : Direction<Dim>::all_directions()) {
377 : fluxes_args_on_faces.emplace(
378 : direction, elliptic::util::apply_at<
379 : domain::make_faces_tags<Dim, fluxes_args_tags,
380 : fluxes_args_volume_tags>,
381 : fluxes_args_volume_tags>(get_items, box, direction));
382 : }
383 : db::mutate<OperatorAppliedToFieldsTag, all_mortar_data_tag>(
384 : [](const auto&... args) {
385 : elliptic::dg::apply_operator<System, Linearized>(args...);
386 : },
387 : make_not_null(&box), db::get<PrimalFieldsTag>(box),
388 : db::get<DerivFieldsTag>(box), db::get<PrimalFluxesTag>(box), element,
389 : db::get<domain::Tags::Mesh<Dim>>(box),
390 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
391 : Frame::Inertial>>(box),
392 : db::get<domain::Tags::DetInvJacobian<Frame::ElementLogical,
393 : Frame::Inertial>>(box),
394 : db::get<
395 : domain::Tags::DetJacobian<Frame::ElementLogical, Frame::Inertial>>(
396 : box),
397 : db::get<domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
398 : Frame::Inertial>>(box),
399 : db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormal<Dim>>>(box),
400 : db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormalVector<Dim>>>(
401 : box),
402 : db::get<domain::Tags::Faces<
403 : Dim, domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>>>(box),
404 : db::get<domain::Tags::Faces<
405 : Dim, domain::Tags::DetSurfaceJacobian<Frame::ElementLogical,
406 : Frame::Inertial>>>(box),
407 : db::get<domain::Tags::Faces<
408 : Dim, domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
409 : Frame::Inertial>>>(box),
410 : db::get<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>>(box),
411 : db::get<::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>>(box),
412 : db::get<::Tags::Mortars<domain::Tags::DetSurfaceJacobian<
413 : Frame::ElementLogical, Frame::Inertial>,
414 : Dim>>(box),
415 : db::get<::Tags::Mortars<elliptic::dg::Tags::PenaltyFactor, Dim>>(box),
416 : db::get<elliptic::dg::Tags::Massive>(box),
417 : db::get<elliptic::dg::Tags::Formulation>(box), temporal_id,
418 : fluxes_args_on_faces,
419 : std::forward_as_tuple(db::get<SourcesArgsTags>(box)...));
420 :
421 : // Increment temporal ID
422 : db::mutate<TemporalIdTag>(
423 : [](const gsl::not_null<size_t*> stored_temporal_id) {
424 : ++(*stored_temporal_id);
425 : },
426 : make_not_null(&box));
427 :
428 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
429 : }
430 : };
431 :
432 : } // namespace detail
433 :
434 : /*!
435 : * \brief Initialize geometric and background quantities for the elliptic DG
436 : * operator
437 : *
438 : * The geometric and background quantities are initialized together because the
439 : * geometry depends on the background metric through the normalization of face
440 : * normals. Other examples for background fields are curvature quantities
441 : * associated with the background metric, or matter sources such as a
442 : * mass-density in the XCTS equations. All `System::background_fields` are
443 : * retrieved from the `BackgroundTag` together, to enable re-using cached
444 : * temporary quantities in the computations. The `variables` function is invoked
445 : * on the `BackgroundTag` with the inertial coordinates, the element's `Mesh`
446 : * and the element's inverse Jacobian. These arguments allow computing numeric
447 : * derivatives, if necessary. The `BackgroundTag` can be set to `void` (default)
448 : * if the `System` has no background fields.
449 : *
450 : * DataBox:
451 : * - Uses:
452 : * - `domain::Tags::InitialExtents<Dim>`
453 : * - `BackgroundTag`
454 : * - Adds:
455 : * - `::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>`
456 : * - `::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>`
457 : * - `::Tags::Variables<background_fields>`
458 : * - Adds on internal and external faces:
459 : * - `domain::Tags::Coordinates<Dim, Frame::Inertial>`
460 : * - `::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<Dim>>`
461 : * - `::Tags::Magnitude<domain::Tags::UnnormalizedFaceNormal<Dim>>`
462 : * - `::Tags::Variables<background_fields>`
463 : *
464 : * \see elliptic::dg::Actions::apply_operator
465 : */
466 : template <typename System, typename BackgroundTag = void>
467 1 : using initialize_operator = tmpl::list<
468 : detail::InitializeFacesMortarsAndBackground<System, BackgroundTag>>;
469 :
470 : /*!
471 : * \brief AMR projectors for the tags added by `initialize_operator`
472 : */
473 : template <typename System, typename BackgroundTag = void>
474 1 : using amr_projectors = tmpl::append<
475 : tmpl::list<elliptic::dg::InitializeFacesAndMortars<
476 : System::volume_dim, typename System::inv_metric_tag, BackgroundTag>>,
477 : tmpl::conditional_t<
478 : std::is_same_v<typename System::background_fields, tmpl::list<>>,
479 : tmpl::list<>,
480 : tmpl::list<elliptic::dg::InitializeBackground<
481 : System::volume_dim, typename System::background_fields,
482 : BackgroundTag>>>>;
483 :
484 : /*!
485 : * \brief Apply the DG operator to the `PrimalFieldsTag` and write the result to
486 : * the `OperatorAppliedToFieldsTag`
487 : *
488 : * Add the `apply_actions` list to the action list of a parallel component to
489 : * compute the elliptic DG operator or its linearization. The operator involves
490 : * a communication between nearest-neighbor elements. See `elliptic::dg` for
491 : * details on the elliptic DG operator. Make sure to add
492 : * `elliptic::dg::Actions::initialize_operator` to the initialization phase of
493 : * your parallel component so the required DataBox tags are set up before
494 : * applying the operator.
495 : *
496 : * The result of the computation is written to the `OperatorAppliedToFieldsTag`.
497 : * Additionally, the primal fluxes are written to the `PrimalFluxesTag` as an
498 : * intermediate result. The auxiliary fields and fluxes are discarded to avoid
499 : * inflating the memory usage.
500 : *
501 : * You can specify the `PrimalMortarFieldsTag` and the `PrimalMortarFluxesTag`
502 : * to re-use mortar-data memory buffers from other operator applications, for
503 : * example when applying the nonlinear and linearized operator. They default to
504 : * the `PrimalFieldsTag` and the `PrimalFluxesTag`, meaning memory buffers
505 : * corresponding to these tags are set up in the DataBox.
506 : *
507 : * \par AMR
508 : * Also add the `amr_projectors` to the list of AMR projectors to support AMR.
509 : */
510 : template <typename System, bool Linearized, typename PrimalFieldsTag,
511 : typename PrimalFluxesTag, typename OperatorAppliedToFieldsTag,
512 : typename PrimalMortarFieldsTag = PrimalFieldsTag,
513 : typename PrimalMortarFluxesTag = PrimalFluxesTag>
514 1 : struct DgOperator {
515 0 : using system = System;
516 0 : using temporal_id_tag = detail::TemporalIdTag;
517 :
518 : private:
519 0 : static constexpr size_t Dim = System::volume_dim;
520 0 : using DerivVarsTag = db::add_tag_prefix<::Tags::deriv, PrimalFieldsTag,
521 : tmpl::size_t<Dim>, Frame::Inertial>;
522 : public:
523 0 : using apply_actions =
524 : tmpl::list<detail::PrepareAndSendMortarData<
525 : System, Linearized, temporal_id_tag, PrimalFieldsTag,
526 : PrimalFluxesTag, OperatorAppliedToFieldsTag,
527 : PrimalMortarFieldsTag, PrimalMortarFluxesTag>,
528 : detail::ReceiveMortarDataAndApplyOperator<
529 : System, Linearized, temporal_id_tag, PrimalFieldsTag,
530 : PrimalFluxesTag, OperatorAppliedToFieldsTag,
531 : PrimalMortarFieldsTag, PrimalMortarFluxesTag>>;
532 0 : using amr_projectors = tmpl::list<
533 : ::amr::projectors::DefaultInitialize<
534 : DerivVarsTag, PrimalFluxesTag, OperatorAppliedToFieldsTag,
535 : ::Tags::Mortars<elliptic::dg::Tags::MortarData<
536 : typename temporal_id_tag::type,
537 : typename PrimalMortarFieldsTag::tags_list,
538 : typename PrimalMortarFluxesTag::tags_list>,
539 : Dim>>,
540 : ::amr::projectors::CopyFromCreatorOrLeaveAsIs<temporal_id_tag>>;
541 : };
542 :
543 : /*!
544 : * \brief For linear systems, impose inhomogeneous boundary conditions as
545 : * contributions to the fixed sources (i.e. the RHS of the equations).
546 : *
547 : * \see elliptic::dg::impose_inhomogeneous_boundary_conditions_on_source
548 : */
549 : template <typename System, typename FixedSourcesTag,
550 : typename FluxesArgsTags =
551 : elliptic::get_fluxes_argument_tags<System, false>,
552 : typename SourcesArgsTags =
553 : elliptic::get_sources_argument_tags<System, false>,
554 : typename ModifyBoundaryDataArgsTags =
555 : elliptic::get_modify_boundary_data_args_tags<System>>
556 1 : struct ImposeInhomogeneousBoundaryConditionsOnSource;
557 :
558 : /// \cond
559 : template <typename System, typename FixedSourcesTag, typename... FluxesArgsTags,
560 : typename... SourcesArgsTags, typename... ModifyBoundaryDataArgsTags>
561 : struct ImposeInhomogeneousBoundaryConditionsOnSource<
562 : System, FixedSourcesTag, tmpl::list<FluxesArgsTags...>,
563 : tmpl::list<SourcesArgsTags...>, tmpl::list<ModifyBoundaryDataArgsTags...>>
564 : : tt::ConformsTo<::amr::protocols::Projector> {
565 : private:
566 : static constexpr size_t Dim = System::volume_dim;
567 : using BoundaryConditionsBase = typename System::boundary_conditions_base;
568 :
569 : public: // DataBox mutator, amr::protocols::Projector
570 : using const_global_cache_tags =
571 : tmpl::list<elliptic::dg::Tags::PenaltyParameter,
572 : elliptic::dg::Tags::Massive, elliptic::dg::Tags::Formulation,
573 : domain::Tags::ExternalBoundaryConditions<Dim>>;
574 :
575 : using return_tags = tmpl::list<::Tags::DataBox>;
576 : using argument_tags = tmpl::list<>;
577 :
578 : template <typename DbTagsList, typename... AmrData>
579 : static void apply(const gsl::not_null<db::DataBox<DbTagsList>*> box_ptr,
580 : const AmrData&... /*amr_data*/) {
581 : // Just to get the same semantics as in actions
582 : auto& box = *box_ptr;
583 : const auto& element_id = db::get<domain::Tags::Element<Dim>>(box).id();
584 : // Used to retrieve items out of the DataBox to forward to functions
585 : const auto get_items = [](const auto&... args) {
586 : return std::forward_as_tuple(args...);
587 : };
588 : const auto& boundary_conditions =
589 : db::get<domain::Tags::ExternalBoundaryConditions<Dim>>(box).at(
590 : element_id.block_id());
591 : const auto apply_boundary_condition = [&box, &boundary_conditions,
592 : &element_id](
593 : const Direction<Dim>& direction,
594 : auto&&... fields_and_fluxes) {
595 : ASSERT(boundary_conditions.contains(direction),
596 : "No boundary condition is available in block "
597 : << element_id.block_id() << " in direction " << direction
598 : << ". Make sure you are setting up boundary conditions when "
599 : "creating the domain.");
600 : ASSERT(dynamic_cast<const BoundaryConditionsBase*>(
601 : boundary_conditions.at(direction).get()) != nullptr,
602 : "The boundary condition in block "
603 : << element_id.block_id() << " in direction " << direction
604 : << " has an unexpected type. Make sure it derives off the "
605 : "'boundary_conditions_base' class set in the system.");
606 : const auto& boundary_condition =
607 : dynamic_cast<const BoundaryConditionsBase&>(
608 : *boundary_conditions.at(direction));
609 : elliptic::apply_boundary_condition<false>(
610 : boundary_condition, box, direction,
611 : std::forward<decltype(fields_and_fluxes)>(fields_and_fluxes)...);
612 : };
613 :
614 : // Can't `db::get` the arguments for the boundary conditions within
615 : // `db::mutate`, so we extract the data to mutate and move it back in when
616 : // we're done.
617 : typename FixedSourcesTag::type fixed_sources;
618 : db::mutate<FixedSourcesTag>(
619 : [&fixed_sources](const auto local_fixed_sources) {
620 : fixed_sources = std::move(*local_fixed_sources);
621 : },
622 : make_not_null(&box));
623 :
624 : using fluxes_args_tags =
625 : typename elliptic::get_fluxes_argument_tags<System, false>;
626 : using fluxes_args_volume_tags =
627 : elliptic::get_fluxes_volume_tags<System, false>;
628 : DirectionMap<Dim, std::tuple<decltype(db::get<FluxesArgsTags>(box))...>>
629 : fluxes_args_on_faces{};
630 : for (const auto& direction : Direction<Dim>::all_directions()) {
631 : fluxes_args_on_faces.emplace(
632 : direction, elliptic::util::apply_at<
633 : domain::make_faces_tags<Dim, fluxes_args_tags,
634 : fluxes_args_volume_tags>,
635 : fluxes_args_volume_tags>(get_items, box, direction));
636 : }
637 :
638 : elliptic::dg::impose_inhomogeneous_boundary_conditions_on_source<System>(
639 : make_not_null(&fixed_sources), db::get<domain::Tags::Element<Dim>>(box),
640 : db::get<domain::Tags::Mesh<Dim>>(box),
641 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
642 : Frame::Inertial>>(box),
643 : db::get<domain::Tags::DetInvJacobian<Frame::ElementLogical,
644 : Frame::Inertial>>(box),
645 : db::get<
646 : domain::Tags::DetJacobian<Frame::ElementLogical, Frame::Inertial>>(
647 : box),
648 : db::get<domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
649 : Frame::Inertial>>(box),
650 : db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormal<Dim>>>(box),
651 : db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormalVector<Dim>>>(
652 : box),
653 : db::get<domain::Tags::Faces<
654 : Dim, domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>>>(box),
655 : db::get<domain::Tags::Faces<
656 : Dim, domain::Tags::DetSurfaceJacobian<Frame::ElementLogical,
657 : Frame::Inertial>>>(box),
658 : db::get<domain::Tags::Faces<
659 : Dim, domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
660 : Frame::Inertial>>>(box),
661 : db::get<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>>(box),
662 : db::get<::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>>(box),
663 : db::get<::Tags::Mortars<domain::Tags::DetSurfaceJacobian<
664 : Frame::ElementLogical, Frame::Inertial>,
665 : Dim>>(box),
666 : db::get<::Tags::Mortars<elliptic::dg::Tags::PenaltyFactor, Dim>>(box),
667 : db::get<elliptic::dg::Tags::Massive>(box),
668 : db::get<elliptic::dg::Tags::Formulation>(box), apply_boundary_condition,
669 : std::forward_as_tuple(db::get<FluxesArgsTags>(box)...),
670 : std::forward_as_tuple(db::get<SourcesArgsTags>(box)...),
671 : fluxes_args_on_faces,
672 : std::forward_as_tuple(db::get<ModifyBoundaryDataArgsTags>(box)...));
673 :
674 : // Move the mutated data back into the DataBox
675 : db::mutate<FixedSourcesTag>(
676 : [&fixed_sources](const auto local_fixed_sources) {
677 : *local_fixed_sources = std::move(fixed_sources);
678 : },
679 : make_not_null(&box));
680 : }
681 : };
682 : /// \endcond
683 :
684 : } // namespace elliptic::dg::Actions
|
