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/TaggedTuple.hpp"
16 : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
17 : #include "Domain/Creators/Tags/ExternalBoundaryConditions.hpp"
18 : #include "Domain/Creators/Tags/InitialExtents.hpp"
19 : #include "Domain/FaceNormal.hpp"
20 : #include "Domain/Structure/Direction.hpp"
21 : #include "Domain/Structure/DirectionMap.hpp"
22 : #include "Domain/Structure/DirectionalIdMap.hpp"
23 : #include "Domain/Structure/Element.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/Initialization.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/Convergence/Tags.hpp"
37 : #include "NumericalAlgorithms/DiscontinuousGalerkin/HasReceivedFromAllMortars.hpp"
38 : #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
39 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags.hpp"
40 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
41 : #include "Parallel/AlgorithmExecution.hpp"
42 : #include "Parallel/GlobalCache.hpp"
43 : #include "Parallel/InboxInserters.hpp"
44 : #include "Parallel/Invoke.hpp"
45 : #include "ParallelAlgorithms/Amr/Projectors/CopyFromCreatorOrLeaveAsIs.hpp"
46 : #include "ParallelAlgorithms/Amr/Projectors/DefaultInitialize.hpp"
47 : #include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
48 : #include "Utilities/ErrorHandling/Assert.hpp"
49 : #include "Utilities/GetOutput.hpp"
50 : #include "Utilities/Gsl.hpp"
51 : #include "Utilities/Literals.hpp"
52 : #include "Utilities/TMPL.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 : typename ModifyBoundaryDataArgsTags =
303 : elliptic::get_modify_boundary_data_args_tags<System, Linearized>>
304 : struct ReceiveMortarDataAndApplyOperator;
305 :
306 : template <typename System, bool Linearized, typename TemporalIdTag,
307 : typename PrimalFieldsTag, typename PrimalFluxesTag,
308 : typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
309 : typename PrimalMortarFluxesTag, typename... FluxesArgsTags,
310 : typename... SourcesArgsTags, typename... ModifyBoundaryDataArgsTags>
311 : struct ReceiveMortarDataAndApplyOperator<
312 : System, Linearized, TemporalIdTag, PrimalFieldsTag, PrimalFluxesTag,
313 : OperatorAppliedToFieldsTag, PrimalMortarFieldsTag, PrimalMortarFluxesTag,
314 : tmpl::list<FluxesArgsTags...>, tmpl::list<SourcesArgsTags...>,
315 : tmpl::list<ModifyBoundaryDataArgsTags...>> {
316 : private:
317 : static constexpr size_t Dim = System::volume_dim;
318 : using all_mortar_data_tag = ::Tags::Mortars<
319 : elliptic::dg::Tags::MortarData<typename TemporalIdTag::type,
320 : typename PrimalMortarFieldsTag::tags_list,
321 : typename PrimalMortarFluxesTag::tags_list>,
322 : Dim>;
323 : using mortar_data_inbox_tag =
324 : MortarDataInboxTag<Dim, TemporalIdTag,
325 : typename PrimalMortarFieldsTag::tags_list,
326 : typename PrimalMortarFluxesTag::tags_list>;
327 : using DerivFieldsTag = db::add_tag_prefix<::Tags::deriv, PrimalFieldsTag,
328 : tmpl::size_t<Dim>, Frame::Inertial>;
329 : public:
330 : using const_global_cache_tags =
331 : tmpl::list<elliptic::dg::Tags::PenaltyParameter,
332 : elliptic::dg::Tags::Massive, elliptic::dg::Tags::Formulation>;
333 : using inbox_tags = tmpl::list<mortar_data_inbox_tag>;
334 :
335 : template <typename DbTags, typename... InboxTags, typename Metavariables,
336 : typename ArrayIndex, typename ActionList,
337 : typename ParallelComponent>
338 : static Parallel::iterable_action_return_t apply(
339 : db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
340 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
341 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
342 : const ParallelComponent* const /*meta*/) {
343 : const auto& temporal_id = get<TemporalIdTag>(box);
344 : const auto& element = get<domain::Tags::Element<Dim>>(box);
345 :
346 : if (not::dg::has_received_from_all_mortars<mortar_data_inbox_tag>(
347 : temporal_id, element, inboxes)) {
348 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
349 : }
350 :
351 : // Move received "remote" mortar data into the DataBox
352 : if (LIKELY(element.number_of_neighbors() > 0)) {
353 : auto received_mortar_data =
354 : std::move(tuples::get<mortar_data_inbox_tag>(inboxes)
355 : .extract(temporal_id)
356 : .mapped());
357 : db::mutate<all_mortar_data_tag>(
358 : [&received_mortar_data, &temporal_id](const auto all_mortar_data) {
359 : for (auto& [mortar_id, mortar_data] : received_mortar_data) {
360 : all_mortar_data->at(mortar_id).remote_insert(
361 : temporal_id, std::move(mortar_data));
362 : }
363 : },
364 : make_not_null(&box));
365 : }
366 :
367 : // Used to retrieve items out of the DataBox to forward to functions
368 : const auto get_items = [](const auto&... args) {
369 : return std::forward_as_tuple(args...);
370 : };
371 :
372 : // Apply DG operator
373 : using fluxes_args_tags =
374 : typename elliptic::get_fluxes_argument_tags<System, Linearized>;
375 : using fluxes_args_volume_tags =
376 : typename elliptic::get_fluxes_volume_tags<System, Linearized>;
377 : DirectionMap<Dim, std::tuple<decltype(db::get<FluxesArgsTags>(box))...>>
378 : fluxes_args_on_faces{};
379 : for (const auto& direction : Direction<Dim>::all_directions()) {
380 : fluxes_args_on_faces.emplace(
381 : direction, elliptic::util::apply_at<
382 : domain::make_faces_tags<Dim, fluxes_args_tags,
383 : fluxes_args_volume_tags>,
384 : fluxes_args_volume_tags>(get_items, box, direction));
385 : }
386 : db::mutate<OperatorAppliedToFieldsTag, all_mortar_data_tag>(
387 : [](const auto&... args) {
388 : elliptic::dg::apply_operator<System, Linearized>(args...);
389 : },
390 : make_not_null(&box), db::get<PrimalFieldsTag>(box),
391 : db::get<DerivFieldsTag>(box), db::get<PrimalFluxesTag>(box), element,
392 : db::get<domain::Tags::Mesh<Dim>>(box),
393 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
394 : Frame::Inertial>>(box),
395 : db::get<domain::Tags::DetInvJacobian<Frame::ElementLogical,
396 : Frame::Inertial>>(box),
397 : db::get<
398 : domain::Tags::DetJacobian<Frame::ElementLogical, Frame::Inertial>>(
399 : box),
400 : db::get<domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
401 : Frame::Inertial>>(box),
402 : db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormal<Dim>>>(box),
403 : db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormalVector<Dim>>>(
404 : box),
405 : db::get<domain::Tags::Faces<
406 : Dim, domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>>>(box),
407 : db::get<domain::Tags::Faces<
408 : Dim, domain::Tags::DetSurfaceJacobian<Frame::ElementLogical,
409 : Frame::Inertial>>>(box),
410 : db::get<domain::Tags::Faces<
411 : Dim, domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
412 : Frame::Inertial>>>(box),
413 : db::get<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>>(box),
414 : db::get<::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>>(box),
415 : db::get<::Tags::Mortars<domain::Tags::DetSurfaceJacobian<
416 : Frame::ElementLogical, Frame::Inertial>,
417 : Dim>>(box),
418 : db::get<::Tags::Mortars<elliptic::dg::Tags::PenaltyFactor, Dim>>(box),
419 : db::get<elliptic::dg::Tags::Massive>(box),
420 : db::get<elliptic::dg::Tags::Formulation>(box), temporal_id,
421 : fluxes_args_on_faces,
422 : std::forward_as_tuple(db::get<SourcesArgsTags>(box)...),
423 : std::forward_as_tuple(db::get<ModifyBoundaryDataArgsTags>(box)...));
424 :
425 : // Increment temporal ID
426 : db::mutate<TemporalIdTag>(
427 : [](const gsl::not_null<size_t*> stored_temporal_id) {
428 : ++(*stored_temporal_id);
429 : },
430 : make_not_null(&box));
431 :
432 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
433 : }
434 : };
435 :
436 : } // namespace detail
437 :
438 : /*!
439 : * \brief Initialize geometric and background quantities for the elliptic DG
440 : * operator
441 : *
442 : * The geometric and background quantities are initialized together because the
443 : * geometry depends on the background metric through the normalization of face
444 : * normals. Other examples for background fields are curvature quantities
445 : * associated with the background metric, or matter sources such as a
446 : * mass-density in the XCTS equations. All `System::background_fields` are
447 : * retrieved from the `BackgroundTag` together, to enable re-using cached
448 : * temporary quantities in the computations. The `variables` function is invoked
449 : * on the `BackgroundTag` with the inertial coordinates, the element's `Mesh`
450 : * and the element's inverse Jacobian. These arguments allow computing numeric
451 : * derivatives, if necessary. The `BackgroundTag` can be set to `void` (default)
452 : * if the `System` has no background fields.
453 : *
454 : * DataBox:
455 : * - Uses:
456 : * - `domain::Tags::InitialExtents<Dim>`
457 : * - `BackgroundTag`
458 : * - Adds:
459 : * - `::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>`
460 : * - `::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>`
461 : * - `::Tags::Variables<background_fields>`
462 : * - Adds on internal and external faces:
463 : * - `domain::Tags::Coordinates<Dim, Frame::Inertial>`
464 : * - `::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<Dim>>`
465 : * - `::Tags::Magnitude<domain::Tags::UnnormalizedFaceNormal<Dim>>`
466 : * - `::Tags::Variables<background_fields>`
467 : *
468 : * \see elliptic::dg::Actions::apply_operator
469 : */
470 : template <typename System, typename BackgroundTag = void>
471 1 : using initialize_operator = tmpl::list<
472 : detail::InitializeFacesMortarsAndBackground<System, BackgroundTag>>;
473 :
474 : /*!
475 : * \brief AMR projectors for the tags added by `initialize_operator`
476 : */
477 : template <typename System, typename BackgroundTag = void>
478 1 : using amr_projectors = tmpl::append<
479 : tmpl::list<elliptic::dg::InitializeFacesAndMortars<
480 : System::volume_dim, typename System::inv_metric_tag, BackgroundTag>>,
481 : tmpl::conditional_t<
482 : std::is_same_v<typename System::background_fields, tmpl::list<>>,
483 : tmpl::list<>,
484 : tmpl::list<elliptic::dg::InitializeBackground<
485 : System::volume_dim, typename System::background_fields,
486 : BackgroundTag>>>>;
487 :
488 : /*!
489 : * \brief Apply the DG operator to the `PrimalFieldsTag` and write the result to
490 : * the `OperatorAppliedToFieldsTag`
491 : *
492 : * Add the `apply_actions` list to the action list of a parallel component to
493 : * compute the elliptic DG operator or its linearization. The operator involves
494 : * a communication between nearest-neighbor elements. See `elliptic::dg` for
495 : * details on the elliptic DG operator. Make sure to add
496 : * `elliptic::dg::Actions::initialize_operator` to the initialization phase of
497 : * your parallel component so the required DataBox tags are set up before
498 : * applying the operator.
499 : *
500 : * The result of the computation is written to the `OperatorAppliedToFieldsTag`.
501 : * Additionally, the primal fluxes are written to the `PrimalFluxesTag` as an
502 : * intermediate result. The auxiliary fields and fluxes are discarded to avoid
503 : * inflating the memory usage.
504 : *
505 : * You can specify the `PrimalMortarFieldsTag` and the `PrimalMortarFluxesTag`
506 : * to re-use mortar-data memory buffers from other operator applications, for
507 : * example when applying the nonlinear and linearized operator. They default to
508 : * the `PrimalFieldsTag` and the `PrimalFluxesTag`, meaning memory buffers
509 : * corresponding to these tags are set up in the DataBox.
510 : *
511 : * \par AMR
512 : * Also add the `amr_projectors` to the list of AMR projectors to support AMR.
513 : */
514 : template <typename System, bool Linearized, typename PrimalFieldsTag,
515 : typename PrimalFluxesTag, typename OperatorAppliedToFieldsTag,
516 : typename PrimalMortarFieldsTag = PrimalFieldsTag,
517 : typename PrimalMortarFluxesTag = PrimalFluxesTag>
518 1 : struct DgOperator {
519 0 : using system = System;
520 0 : using temporal_id_tag = detail::TemporalIdTag;
521 :
522 : private:
523 0 : static constexpr size_t Dim = System::volume_dim;
524 0 : using DerivVarsTag = db::add_tag_prefix<::Tags::deriv, PrimalFieldsTag,
525 : tmpl::size_t<Dim>, Frame::Inertial>;
526 : public:
527 0 : using apply_actions =
528 : tmpl::list<detail::PrepareAndSendMortarData<
529 : System, Linearized, temporal_id_tag, PrimalFieldsTag,
530 : PrimalFluxesTag, OperatorAppliedToFieldsTag,
531 : PrimalMortarFieldsTag, PrimalMortarFluxesTag>,
532 : detail::ReceiveMortarDataAndApplyOperator<
533 : System, Linearized, temporal_id_tag, PrimalFieldsTag,
534 : PrimalFluxesTag, OperatorAppliedToFieldsTag,
535 : PrimalMortarFieldsTag, PrimalMortarFluxesTag>>;
536 0 : using amr_projectors = tmpl::list<
537 : ::amr::projectors::DefaultInitialize<
538 : DerivVarsTag, PrimalFluxesTag, OperatorAppliedToFieldsTag,
539 : ::Tags::Mortars<elliptic::dg::Tags::MortarData<
540 : typename temporal_id_tag::type,
541 : typename PrimalMortarFieldsTag::tags_list,
542 : typename PrimalMortarFluxesTag::tags_list>,
543 : Dim>>,
544 : ::amr::projectors::CopyFromCreatorOrLeaveAsIs<temporal_id_tag>>;
545 : };
546 :
547 : /*!
548 : * \brief For linear systems, impose inhomogeneous boundary conditions as
549 : * contributions to the fixed sources (i.e. the RHS of the equations).
550 : *
551 : * \see elliptic::dg::impose_inhomogeneous_boundary_conditions_on_source
552 : */
553 : template <typename System, typename FixedSourcesTag,
554 : typename FluxesArgsTags =
555 : elliptic::get_fluxes_argument_tags<System, false>,
556 : typename SourcesArgsTags =
557 : elliptic::get_sources_argument_tags<System, false>,
558 : typename ModifyBoundaryDataArgsTags =
559 : elliptic::get_modify_boundary_data_args_tags<System, false>>
560 1 : struct ImposeInhomogeneousBoundaryConditionsOnSource;
561 :
562 : /// \cond
563 : template <typename System, typename FixedSourcesTag, typename... FluxesArgsTags,
564 : typename... SourcesArgsTags, typename... ModifyBoundaryDataArgsTags>
565 : struct ImposeInhomogeneousBoundaryConditionsOnSource<
566 : System, FixedSourcesTag, tmpl::list<FluxesArgsTags...>,
567 : tmpl::list<SourcesArgsTags...>, tmpl::list<ModifyBoundaryDataArgsTags...>>
568 : : tt::ConformsTo<::amr::protocols::Projector> {
569 : private:
570 : static constexpr size_t Dim = System::volume_dim;
571 : using BoundaryConditionsBase = typename System::boundary_conditions_base;
572 :
573 : public: // DataBox mutator, amr::protocols::Projector
574 : using const_global_cache_tags =
575 : tmpl::list<elliptic::dg::Tags::PenaltyParameter,
576 : elliptic::dg::Tags::Massive, elliptic::dg::Tags::Formulation,
577 : domain::Tags::ExternalBoundaryConditions<Dim>>;
578 :
579 : using return_tags = tmpl::list<::Tags::DataBox>;
580 : using argument_tags = tmpl::list<>;
581 :
582 : template <typename DbTagsList, typename... AmrData>
583 : static void apply(const gsl::not_null<db::DataBox<DbTagsList>*> box_ptr,
584 : const AmrData&... /*amr_data*/) {
585 : // Just to get the same semantics as in actions
586 : auto& box = *box_ptr;
587 : const auto& element_id = db::get<domain::Tags::Element<Dim>>(box).id();
588 : // Used to retrieve items out of the DataBox to forward to functions
589 : const auto get_items = [](const auto&... args) {
590 : return std::forward_as_tuple(args...);
591 : };
592 : const auto& boundary_conditions =
593 : db::get<domain::Tags::ExternalBoundaryConditions<Dim>>(box).at(
594 : element_id.block_id());
595 : const auto apply_boundary_condition = [&box, &boundary_conditions,
596 : &element_id](
597 : const Direction<Dim>& direction,
598 : auto&&... fields_and_fluxes) {
599 : ASSERT(boundary_conditions.contains(direction),
600 : "No boundary condition is available in block "
601 : << element_id.block_id() << " in direction " << direction
602 : << ". Make sure you are setting up boundary conditions when "
603 : "creating the domain.");
604 : ASSERT(dynamic_cast<const BoundaryConditionsBase*>(
605 : boundary_conditions.at(direction).get()) != nullptr,
606 : "The boundary condition in block "
607 : << element_id.block_id() << " in direction " << direction
608 : << " has an unexpected type. Make sure it derives off the "
609 : "'boundary_conditions_base' class set in the system.");
610 : const auto& boundary_condition =
611 : dynamic_cast<const BoundaryConditionsBase&>(
612 : *boundary_conditions.at(direction));
613 : elliptic::apply_boundary_condition<false>(
614 : boundary_condition, box, direction,
615 : std::forward<decltype(fields_and_fluxes)>(fields_and_fluxes)...);
616 : };
617 :
618 : // Can't `db::get` the arguments for the boundary conditions within
619 : // `db::mutate`, so we extract the data to mutate and move it back in when
620 : // we're done.
621 : typename FixedSourcesTag::type fixed_sources;
622 : db::mutate<FixedSourcesTag>(
623 : [&fixed_sources](const auto local_fixed_sources) {
624 : fixed_sources = std::move(*local_fixed_sources);
625 : },
626 : make_not_null(&box));
627 :
628 : using fluxes_args_tags =
629 : typename elliptic::get_fluxes_argument_tags<System, false>;
630 : using fluxes_args_volume_tags =
631 : elliptic::get_fluxes_volume_tags<System, false>;
632 : DirectionMap<Dim, std::tuple<decltype(db::get<FluxesArgsTags>(box))...>>
633 : fluxes_args_on_faces{};
634 : for (const auto& direction : Direction<Dim>::all_directions()) {
635 : fluxes_args_on_faces.emplace(
636 : direction, elliptic::util::apply_at<
637 : domain::make_faces_tags<Dim, fluxes_args_tags,
638 : fluxes_args_volume_tags>,
639 : fluxes_args_volume_tags>(get_items, box, direction));
640 : }
641 :
642 : elliptic::dg::impose_inhomogeneous_boundary_conditions_on_source<System>(
643 : make_not_null(&fixed_sources), db::get<domain::Tags::Element<Dim>>(box),
644 : db::get<domain::Tags::Mesh<Dim>>(box),
645 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
646 : Frame::Inertial>>(box),
647 : db::get<domain::Tags::DetInvJacobian<Frame::ElementLogical,
648 : Frame::Inertial>>(box),
649 : db::get<
650 : domain::Tags::DetJacobian<Frame::ElementLogical, Frame::Inertial>>(
651 : box),
652 : db::get<domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
653 : Frame::Inertial>>(box),
654 : db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormal<Dim>>>(box),
655 : db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormalVector<Dim>>>(
656 : box),
657 : db::get<domain::Tags::Faces<
658 : Dim, domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>>>(box),
659 : db::get<domain::Tags::Faces<
660 : Dim, domain::Tags::DetSurfaceJacobian<Frame::ElementLogical,
661 : Frame::Inertial>>>(box),
662 : db::get<domain::Tags::Faces<
663 : Dim, domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
664 : Frame::Inertial>>>(box),
665 : db::get<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>>(box),
666 : db::get<::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>>(box),
667 : db::get<::Tags::Mortars<domain::Tags::DetSurfaceJacobian<
668 : Frame::ElementLogical, Frame::Inertial>,
669 : Dim>>(box),
670 : db::get<::Tags::Mortars<elliptic::dg::Tags::PenaltyFactor, Dim>>(box),
671 : db::get<elliptic::dg::Tags::Massive>(box),
672 : db::get<elliptic::dg::Tags::Formulation>(box), apply_boundary_condition,
673 : std::forward_as_tuple(db::get<FluxesArgsTags>(box)...),
674 : std::forward_as_tuple(db::get<SourcesArgsTags>(box)...),
675 : fluxes_args_on_faces,
676 : std::forward_as_tuple(db::get<ModifyBoundaryDataArgsTags>(box)...));
677 :
678 : // Move the mutated data back into the DataBox
679 : db::mutate<FixedSourcesTag>(
680 : [&fixed_sources](const auto local_fixed_sources) {
681 : *local_fixed_sources = std::move(fixed_sources);
682 : },
683 : make_not_null(&box));
684 : }
685 : };
686 : /// \endcond
687 :
688 : } // namespace elliptic::dg::Actions
|
