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 :
164 : public:
165 : // Request these tags be added to the DataBox. We
166 : // don't actually need to initialize them, because the
167 : // `PrimalFieldsTag` will be set by other actions before applying the operator
168 : // and the remaining tags hold output of the operator.
169 : using simple_tags =
170 : tmpl::list<TemporalIdTag, PrimalFieldsTag, PrimalFluxesTag,
171 : OperatorAppliedToFieldsTag, all_mortar_data_tag>;
172 : using compute_tags = tmpl::list<>;
173 : using const_global_cache_tags =
174 : tmpl::list<domain::Tags::ExternalBoundaryConditions<Dim>>;
175 :
176 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
177 : typename ActionList, typename ParallelComponent>
178 : static Parallel::iterable_action_return_t apply(
179 : db::DataBox<DbTagsList>& box,
180 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
181 : Parallel::GlobalCache<Metavariables>& cache,
182 : const ElementId<Dim>& element_id, const ActionList /*meta*/,
183 : const ParallelComponent* const /*meta*/) {
184 : const auto& temporal_id = db::get<TemporalIdTag>(box);
185 : const auto& element = db::get<domain::Tags::Element<Dim>>(box);
186 : const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
187 : const size_t num_points = mesh.number_of_grid_points();
188 : const auto& mortar_meshes =
189 : db::get<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>>(box);
190 : const auto& boundary_conditions =
191 : db::get<domain::Tags::ExternalBoundaryConditions<Dim>>(box).at(
192 : element_id.block_id());
193 : const auto apply_boundary_condition = [&box, &boundary_conditions,
194 : &element_id](
195 : const Direction<Dim>& direction,
196 : auto&&... fields_and_fluxes) {
197 : ASSERT(boundary_conditions.contains(direction),
198 : "No boundary condition is available in block "
199 : << element_id.block_id() << " in direction " << direction
200 : << ". Make sure you are setting up boundary conditions when "
201 : "creating the domain.");
202 : ASSERT(dynamic_cast<const BoundaryConditionsBase*>(
203 : boundary_conditions.at(direction).get()) != nullptr,
204 : "The boundary condition in block "
205 : << element_id.block_id() << " in direction " << direction
206 : << " has an unexpected type. Make sure it derives off the "
207 : "'boundary_conditions_base' class set in the system.");
208 : const auto& boundary_condition =
209 : dynamic_cast<const BoundaryConditionsBase&>(
210 : *boundary_conditions.at(direction));
211 : elliptic::apply_boundary_condition<Linearized>(
212 : boundary_condition, box, direction,
213 : std::forward<decltype(fields_and_fluxes)>(fields_and_fluxes)...);
214 : };
215 :
216 : // Can't `db::get` the arguments for the boundary conditions within
217 : // `db::mutate`, so we extract the data to mutate and move it back in when
218 : // we're done.
219 : // Possible optimization: also keep memory for the mortar data around in the
220 : // DataBox. Currently the mortar data is extracted and erased by
221 : // `apply_operator` anyway, so we can just create a new map here to avoid
222 : // dealing with AMR resizing the mortars. When we keep the memory around,
223 : // its size has to be adjusted when the mesh changes during AMR.
224 : typename PrimalFluxesTag::type primal_fluxes;
225 : typename all_mortar_data_tag::type all_mortar_data{};
226 : db::mutate<PrimalFluxesTag>(
227 : [&primal_fluxes](const auto local_primal_fluxes) {
228 : primal_fluxes = std::move(*local_primal_fluxes);
229 : },
230 : make_not_null(&box));
231 :
232 : // Prepare mortar data
233 : //
234 : // These memory buffers will be discarded when the action returns so we
235 : // don't inflate the memory usage of the simulation when the element is
236 : // inactive.
237 : Variables<db::wrap_tags_in<::Tags::deriv,
238 : typename PrimalFieldsTag::type::tags_list,
239 : tmpl::size_t<Dim>, Frame::Inertial>>
240 : deriv_fields{num_points};
241 : elliptic::dg::prepare_mortar_data<System, Linearized>(
242 : make_not_null(&deriv_fields), make_not_null(&primal_fluxes),
243 : make_not_null(&all_mortar_data), db::get<PrimalFieldsTag>(box), element,
244 : db::get<domain::Tags::Mesh<Dim>>(box),
245 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
246 : Frame::Inertial>>(box),
247 : db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormal<Dim>>>(box),
248 : mortar_meshes,
249 : db::get<::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>>(box),
250 : temporal_id, apply_boundary_condition,
251 : std::forward_as_tuple(db::get<FluxesArgsTags>(box)...));
252 :
253 : // Move the mutated data back into the DataBox
254 : db::mutate<PrimalFluxesTag, all_mortar_data_tag>(
255 : [&primal_fluxes, &all_mortar_data](const auto local_primal_fluxes,
256 : const auto local_all_mortar_data) {
257 : *local_primal_fluxes = std::move(primal_fluxes);
258 : *local_all_mortar_data = std::move(all_mortar_data);
259 : },
260 : make_not_null(&box));
261 :
262 : // Send mortar data to neighbors
263 : auto& receiver_proxy =
264 : Parallel::get_parallel_component<ParallelComponent>(cache);
265 : for (const auto& [direction, neighbors] : element.neighbors()) {
266 : const size_t dimension = direction.dimension();
267 : for (const auto& neighbor_id : neighbors) {
268 : const auto& orientation = neighbors.orientation(neighbor_id);
269 : const auto direction_from_neighbor = orientation(direction.opposite());
270 : const ::dg::MortarId<Dim> mortar_id{direction, neighbor_id};
271 : // Make a copy of the local boundary data on the mortar to send to the
272 : // neighbor
273 : auto remote_boundary_data_on_mortar =
274 : get<all_mortar_data_tag>(box).at(mortar_id).local_data(temporal_id);
275 : // Reorient the data to the neighbor orientation if necessary
276 : if (not orientation.is_aligned()) {
277 : remote_boundary_data_on_mortar.orient_on_slice(
278 : mortar_meshes.at(mortar_id).extents(), dimension, orientation);
279 : }
280 : // Send remote data to neighbor
281 : Parallel::receive_data<mortar_data_inbox_tag>(
282 : receiver_proxy[neighbor_id], temporal_id,
283 : std::make_pair(
284 : ::dg::MortarId<Dim>{direction_from_neighbor, element.id()},
285 : std::move(remote_boundary_data_on_mortar)));
286 : } // loop over neighbors in direction
287 : } // loop over directions
288 :
289 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
290 : }
291 : };
292 :
293 : // Wait until all mortar data from neighbors is available. Then add boundary
294 : // corrections to the primal fluxes, compute their derivatives (i.e. the second
295 : // derivatives of the primal variables) and add boundary corrections to the
296 : // result.
297 : template <typename System, bool Linearized, typename TemporalIdTag,
298 : typename PrimalFieldsTag, typename PrimalFluxesTag,
299 : typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
300 : typename PrimalMortarFluxesTag,
301 : typename FluxesArgsTags =
302 : elliptic::get_fluxes_argument_tags<System, Linearized>,
303 : typename SourcesArgsTags =
304 : elliptic::get_sources_argument_tags<System, Linearized>>
305 : struct ReceiveMortarDataAndApplyOperator;
306 :
307 : template <typename System, bool Linearized, typename TemporalIdTag,
308 : typename PrimalFieldsTag, typename PrimalFluxesTag,
309 : typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
310 : typename PrimalMortarFluxesTag, typename... FluxesArgsTags,
311 : typename... SourcesArgsTags>
312 : struct ReceiveMortarDataAndApplyOperator<
313 : System, Linearized, TemporalIdTag, PrimalFieldsTag, PrimalFluxesTag,
314 : OperatorAppliedToFieldsTag, PrimalMortarFieldsTag, PrimalMortarFluxesTag,
315 : tmpl::list<FluxesArgsTags...>, tmpl::list<SourcesArgsTags...>> {
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 :
328 : public:
329 : using const_global_cache_tags =
330 : tmpl::list<elliptic::dg::Tags::PenaltyParameter,
331 : elliptic::dg::Tags::Massive, elliptic::dg::Tags::Formulation>;
332 : using inbox_tags = tmpl::list<mortar_data_inbox_tag>;
333 :
334 : template <typename DbTags, typename... InboxTags, typename Metavariables,
335 : typename ArrayIndex, typename ActionList,
336 : typename ParallelComponent>
337 : static Parallel::iterable_action_return_t apply(
338 : db::DataBox<DbTags>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
339 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
340 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
341 : const ParallelComponent* const /*meta*/) {
342 : const auto& temporal_id = get<TemporalIdTag>(box);
343 : const auto& element = get<domain::Tags::Element<Dim>>(box);
344 :
345 : if (not ::dg::has_received_from_all_mortars<mortar_data_inbox_tag>(
346 : temporal_id, element, inboxes)) {
347 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
348 : }
349 :
350 : // Move received "remote" mortar data into the DataBox
351 : if (LIKELY(element.number_of_neighbors() > 0)) {
352 : auto received_mortar_data =
353 : std::move(tuples::get<mortar_data_inbox_tag>(inboxes)
354 : .extract(temporal_id)
355 : .mapped());
356 : db::mutate<all_mortar_data_tag>(
357 : [&received_mortar_data, &temporal_id](const auto all_mortar_data) {
358 : for (auto& [mortar_id, mortar_data] : received_mortar_data) {
359 : all_mortar_data->at(mortar_id).remote_insert(
360 : temporal_id, std::move(mortar_data));
361 : }
362 : },
363 : make_not_null(&box));
364 : }
365 :
366 : // Used to retrieve items out of the DataBox to forward to functions
367 : const auto get_items = [](const auto&... args) {
368 : return std::forward_as_tuple(args...);
369 : };
370 :
371 : // Apply DG operator
372 : using fluxes_args_tags =
373 : typename elliptic::get_fluxes_argument_tags<System, Linearized>;
374 : using fluxes_args_volume_tags =
375 : typename elliptic::get_fluxes_volume_tags<System, Linearized>;
376 : DirectionMap<Dim, std::tuple<decltype(db::get<FluxesArgsTags>(box))...>>
377 : fluxes_args_on_faces{};
378 : for (const auto& direction : Direction<Dim>::all_directions()) {
379 : fluxes_args_on_faces.emplace(
380 : direction, elliptic::util::apply_at<
381 : domain::make_faces_tags<Dim, fluxes_args_tags,
382 : fluxes_args_volume_tags>,
383 : fluxes_args_volume_tags>(get_items, box, direction));
384 : }
385 : db::mutate<OperatorAppliedToFieldsTag, all_mortar_data_tag>(
386 : [](const auto&... args) {
387 : elliptic::dg::apply_operator<System, Linearized>(args...);
388 : },
389 : make_not_null(&box), db::get<PrimalFieldsTag>(box),
390 : db::get<PrimalFluxesTag>(box), element,
391 : db::get<domain::Tags::Mesh<Dim>>(box),
392 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
393 : Frame::Inertial>>(box),
394 : db::get<domain::Tags::DetInvJacobian<Frame::ElementLogical,
395 : Frame::Inertial>>(box),
396 : db::get<
397 : domain::Tags::DetJacobian<Frame::ElementLogical, Frame::Inertial>>(
398 : box),
399 : db::get<domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
400 : Frame::Inertial>>(box),
401 : db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormal<Dim>>>(box),
402 : db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormalVector<Dim>>>(
403 : box),
404 : db::get<domain::Tags::Faces<
405 : Dim, domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>>>(box),
406 : db::get<domain::Tags::Faces<
407 : Dim, domain::Tags::DetSurfaceJacobian<Frame::ElementLogical,
408 : Frame::Inertial>>>(box),
409 : db::get<domain::Tags::Faces<
410 : Dim, domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
411 : Frame::Inertial>>>(box),
412 : db::get<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>>(box),
413 : db::get<::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>>(box),
414 : db::get<::Tags::Mortars<domain::Tags::DetSurfaceJacobian<
415 : Frame::ElementLogical, Frame::Inertial>,
416 : Dim>>(box),
417 : db::get<::Tags::Mortars<elliptic::dg::Tags::PenaltyFactor, Dim>>(box),
418 : db::get<elliptic::dg::Tags::Massive>(box),
419 : db::get<elliptic::dg::Tags::Formulation>(box), temporal_id,
420 : fluxes_args_on_faces,
421 : std::forward_as_tuple(db::get<SourcesArgsTags>(box)...));
422 :
423 : // Increment temporal ID
424 : db::mutate<TemporalIdTag>(
425 : [](const gsl::not_null<size_t*> stored_temporal_id) {
426 : ++(*stored_temporal_id);
427 : },
428 : make_not_null(&box));
429 :
430 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
431 : }
432 : };
433 :
434 : } // namespace detail
435 :
436 : /*!
437 : * \brief Initialize geometric and background quantities for the elliptic DG
438 : * operator
439 : *
440 : * The geometric and background quantities are initialized together because the
441 : * geometry depends on the background metric through the normalization of face
442 : * normals. Other examples for background fields are curvature quantities
443 : * associated with the background metric, or matter sources such as a
444 : * mass-density in the XCTS equations. All `System::background_fields` are
445 : * retrieved from the `BackgroundTag` together, to enable re-using cached
446 : * temporary quantities in the computations. The `variables` function is invoked
447 : * on the `BackgroundTag` with the inertial coordinates, the element's `Mesh`
448 : * and the element's inverse Jacobian. These arguments allow computing numeric
449 : * derivatives, if necessary. The `BackgroundTag` can be set to `void` (default)
450 : * if the `System` has no background fields.
451 : *
452 : * DataBox:
453 : * - Uses:
454 : * - `domain::Tags::InitialExtents<Dim>`
455 : * - `BackgroundTag`
456 : * - Adds:
457 : * - `::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>`
458 : * - `::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>`
459 : * - `::Tags::Variables<background_fields>`
460 : * - Adds on internal and external faces:
461 : * - `domain::Tags::Coordinates<Dim, Frame::Inertial>`
462 : * - `::Tags::Normalized<domain::Tags::UnnormalizedFaceNormal<Dim>>`
463 : * - `::Tags::Magnitude<domain::Tags::UnnormalizedFaceNormal<Dim>>`
464 : * - `::Tags::Variables<background_fields>`
465 : *
466 : * \see elliptic::dg::Actions::apply_operator
467 : */
468 : template <typename System, typename BackgroundTag = void>
469 1 : using initialize_operator = tmpl::list<
470 : detail::InitializeFacesMortarsAndBackground<System, BackgroundTag>>;
471 :
472 : /*!
473 : * \brief AMR projectors for the tags added by `initialize_operator`
474 : */
475 : template <typename System, typename BackgroundTag = void>
476 1 : using amr_projectors = tmpl::append<
477 : tmpl::list<elliptic::dg::InitializeFacesAndMortars<
478 : System::volume_dim, typename System::inv_metric_tag, BackgroundTag>>,
479 : tmpl::conditional_t<
480 : std::is_same_v<typename System::background_fields, tmpl::list<>>,
481 : tmpl::list<>,
482 : tmpl::list<elliptic::dg::InitializeBackground<
483 : System::volume_dim, typename System::background_fields,
484 : BackgroundTag>>>>;
485 :
486 : /*!
487 : * \brief Apply the DG operator to the `PrimalFieldsTag` and write the result to
488 : * the `OperatorAppliedToFieldsTag`
489 : *
490 : * Add the `apply_actions` list to the action list of a parallel component to
491 : * compute the elliptic DG operator or its linearization. The operator involves
492 : * a communication between nearest-neighbor elements. See `elliptic::dg` for
493 : * details on the elliptic DG operator. Make sure to add
494 : * `elliptic::dg::Actions::initialize_operator` to the initialization phase of
495 : * your parallel component so the required DataBox tags are set up before
496 : * applying the operator.
497 : *
498 : * The result of the computation is written to the `OperatorAppliedToFieldsTag`.
499 : * Additionally, the primal fluxes are written to the `PrimalFluxesTag` as an
500 : * intermediate result. The auxiliary fields and fluxes are discarded to avoid
501 : * inflating the memory usage.
502 : *
503 : * You can specify the `PrimalMortarFieldsTag` and the `PrimalMortarFluxesTag`
504 : * to re-use mortar-data memory buffers from other operator applications, for
505 : * example when applying the nonlinear and linearized operator. They default to
506 : * the `PrimalFieldsTag` and the `PrimalFluxesTag`, meaning memory buffers
507 : * corresponding to these tags are set up in the DataBox.
508 : *
509 : * \par AMR
510 : * Also add the `amr_projectors` to the list of AMR projectors to support AMR.
511 : */
512 : template <typename System, bool Linearized, typename PrimalFieldsTag,
513 : typename PrimalFluxesTag, typename OperatorAppliedToFieldsTag,
514 : typename PrimalMortarFieldsTag = PrimalFieldsTag,
515 : typename PrimalMortarFluxesTag = PrimalFluxesTag>
516 1 : struct DgOperator {
517 0 : using system = System;
518 0 : using temporal_id_tag = detail::TemporalIdTag;
519 :
520 : private:
521 0 : static constexpr size_t Dim = System::volume_dim;
522 :
523 : public:
524 0 : using apply_actions =
525 : tmpl::list<detail::PrepareAndSendMortarData<
526 : System, Linearized, temporal_id_tag, PrimalFieldsTag,
527 : PrimalFluxesTag, OperatorAppliedToFieldsTag,
528 : PrimalMortarFieldsTag, PrimalMortarFluxesTag>,
529 : detail::ReceiveMortarDataAndApplyOperator<
530 : System, Linearized, temporal_id_tag, PrimalFieldsTag,
531 : PrimalFluxesTag, OperatorAppliedToFieldsTag,
532 : PrimalMortarFieldsTag, PrimalMortarFluxesTag>>;
533 0 : using amr_projectors = tmpl::list<
534 : ::amr::projectors::DefaultInitialize<
535 : PrimalFluxesTag, OperatorAppliedToFieldsTag,
536 : ::Tags::Mortars<elliptic::dg::Tags::MortarData<
537 : typename temporal_id_tag::type,
538 : typename PrimalMortarFieldsTag::tags_list,
539 : typename PrimalMortarFluxesTag::tags_list>,
540 : Dim>>,
541 : ::amr::projectors::CopyFromCreatorOrLeaveAsIs<temporal_id_tag>>;
542 : };
543 :
544 : /*!
545 : * \brief For linear systems, impose inhomogeneous boundary conditions as
546 : * contributions to the fixed sources (i.e. the RHS of the equations).
547 : *
548 : * \see elliptic::dg::impose_inhomogeneous_boundary_conditions_on_source
549 : */
550 : template <typename System, typename FixedSourcesTag,
551 : typename FluxesArgsTags =
552 : elliptic::get_fluxes_argument_tags<System, false>,
553 : typename SourcesArgsTags =
554 : elliptic::get_sources_argument_tags<System, false>,
555 : typename ModifyBoundaryDataArgsTags =
556 : elliptic::get_modify_boundary_data_args_tags<System>>
557 1 : struct ImposeInhomogeneousBoundaryConditionsOnSource;
558 :
559 : /// \cond
560 : template <typename System, typename FixedSourcesTag, typename... FluxesArgsTags,
561 : typename... SourcesArgsTags, typename... ModifyBoundaryDataArgsTags>
562 : struct ImposeInhomogeneousBoundaryConditionsOnSource<
563 : System, FixedSourcesTag, tmpl::list<FluxesArgsTags...>,
564 : tmpl::list<SourcesArgsTags...>, tmpl::list<ModifyBoundaryDataArgsTags...>>
565 : : tt::ConformsTo<::amr::protocols::Projector> {
566 : private:
567 : static constexpr size_t Dim = System::volume_dim;
568 : using BoundaryConditionsBase = typename System::boundary_conditions_base;
569 :
570 : public: // DataBox mutator, amr::protocols::Projector
571 : using const_global_cache_tags =
572 : tmpl::list<elliptic::dg::Tags::PenaltyParameter,
573 : elliptic::dg::Tags::Massive, elliptic::dg::Tags::Formulation,
574 : domain::Tags::ExternalBoundaryConditions<Dim>>;
575 :
576 : using return_tags = tmpl::list<::Tags::DataBox>;
577 : using argument_tags = tmpl::list<>;
578 :
579 : template <typename DbTagsList, typename... AmrData>
580 : static void apply(const gsl::not_null<db::DataBox<DbTagsList>*> box_ptr,
581 : const AmrData&... /*amr_data*/) {
582 : // Just to get the same semantics as in actions
583 : auto& box = *box_ptr;
584 : const auto& element_id = db::get<domain::Tags::Element<Dim>>(box).id();
585 : // Used to retrieve items out of the DataBox to forward to functions
586 : const auto get_items = [](const auto&... args) {
587 : return std::forward_as_tuple(args...);
588 : };
589 : const auto& boundary_conditions =
590 : db::get<domain::Tags::ExternalBoundaryConditions<Dim>>(box).at(
591 : element_id.block_id());
592 : const auto apply_boundary_condition = [&box, &boundary_conditions,
593 : &element_id](
594 : const Direction<Dim>& direction,
595 : auto&&... fields_and_fluxes) {
596 : ASSERT(boundary_conditions.contains(direction),
597 : "No boundary condition is available in block "
598 : << element_id.block_id() << " in direction " << direction
599 : << ". Make sure you are setting up boundary conditions when "
600 : "creating the domain.");
601 : ASSERT(dynamic_cast<const BoundaryConditionsBase*>(
602 : boundary_conditions.at(direction).get()) != nullptr,
603 : "The boundary condition in block "
604 : << element_id.block_id() << " in direction " << direction
605 : << " has an unexpected type. Make sure it derives off the "
606 : "'boundary_conditions_base' class set in the system.");
607 : const auto& boundary_condition =
608 : dynamic_cast<const BoundaryConditionsBase&>(
609 : *boundary_conditions.at(direction));
610 : elliptic::apply_boundary_condition<false>(
611 : boundary_condition, box, direction,
612 : std::forward<decltype(fields_and_fluxes)>(fields_and_fluxes)...);
613 : };
614 :
615 : // Can't `db::get` the arguments for the boundary conditions within
616 : // `db::mutate`, so we extract the data to mutate and move it back in when
617 : // we're done.
618 : typename FixedSourcesTag::type fixed_sources;
619 : db::mutate<FixedSourcesTag>(
620 : [&fixed_sources](const auto local_fixed_sources) {
621 : fixed_sources = std::move(*local_fixed_sources);
622 : },
623 : make_not_null(&box));
624 :
625 : using fluxes_args_tags =
626 : typename elliptic::get_fluxes_argument_tags<System, false>;
627 : using fluxes_args_volume_tags =
628 : elliptic::get_fluxes_volume_tags<System, false>;
629 : DirectionMap<Dim, std::tuple<decltype(db::get<FluxesArgsTags>(box))...>>
630 : fluxes_args_on_faces{};
631 : for (const auto& direction : Direction<Dim>::all_directions()) {
632 : fluxes_args_on_faces.emplace(
633 : direction, elliptic::util::apply_at<
634 : domain::make_faces_tags<Dim, fluxes_args_tags,
635 : fluxes_args_volume_tags>,
636 : fluxes_args_volume_tags>(get_items, box, direction));
637 : }
638 :
639 : elliptic::dg::impose_inhomogeneous_boundary_conditions_on_source<System>(
640 : make_not_null(&fixed_sources), db::get<domain::Tags::Element<Dim>>(box),
641 : db::get<domain::Tags::Mesh<Dim>>(box),
642 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
643 : Frame::Inertial>>(box),
644 : db::get<domain::Tags::DetInvJacobian<Frame::ElementLogical,
645 : Frame::Inertial>>(box),
646 : db::get<
647 : domain::Tags::DetJacobian<Frame::ElementLogical, Frame::Inertial>>(
648 : box),
649 : db::get<domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
650 : Frame::Inertial>>(box),
651 : db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormal<Dim>>>(box),
652 : db::get<domain::Tags::Faces<Dim, domain::Tags::FaceNormalVector<Dim>>>(
653 : box),
654 : db::get<domain::Tags::Faces<
655 : Dim, domain::Tags::UnnormalizedFaceNormalMagnitude<Dim>>>(box),
656 : db::get<domain::Tags::Faces<
657 : Dim, domain::Tags::DetSurfaceJacobian<Frame::ElementLogical,
658 : Frame::Inertial>>>(box),
659 : db::get<domain::Tags::Faces<
660 : Dim, domain::Tags::DetTimesInvJacobian<Dim, Frame::ElementLogical,
661 : Frame::Inertial>>>(box),
662 : db::get<::Tags::Mortars<domain::Tags::Mesh<Dim - 1>, Dim>>(box),
663 : db::get<::Tags::Mortars<::Tags::MortarSize<Dim - 1>, Dim>>(box),
664 : db::get<::Tags::Mortars<domain::Tags::DetSurfaceJacobian<
665 : Frame::ElementLogical, Frame::Inertial>,
666 : Dim>>(box),
667 : db::get<::Tags::Mortars<elliptic::dg::Tags::PenaltyFactor, Dim>>(box),
668 : db::get<elliptic::dg::Tags::Massive>(box),
669 : db::get<elliptic::dg::Tags::Formulation>(box), apply_boundary_condition,
670 : std::forward_as_tuple(db::get<FluxesArgsTags>(box)...),
671 : std::forward_as_tuple(db::get<SourcesArgsTags>(box)...),
672 : fluxes_args_on_faces,
673 : std::forward_as_tuple(db::get<ModifyBoundaryDataArgsTags>(box)...));
674 :
675 : // Move the mutated data back into the DataBox
676 : db::mutate<FixedSourcesTag>(
677 : [&fixed_sources](const auto local_fixed_sources) {
678 : *local_fixed_sources = std::move(fixed_sources);
679 : },
680 : make_not_null(&box));
681 : }
682 : };
683 : /// \endcond
684 :
685 : } // namespace elliptic::dg::Actions
|
