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