15 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
20 #include "Evolution/BoundaryCorrectionTags.hpp"
21 #include "Evolution/DiscontinuousGalerkin/InboxTags.hpp"
22 #include "Evolution/DiscontinuousGalerkin/LiftFromBoundary.hpp"
23 #include "Evolution/DiscontinuousGalerkin/MortarData.hpp"
24 #include "Evolution/DiscontinuousGalerkin/MortarTags.hpp"
25 #include "Evolution/DiscontinuousGalerkin/NormalVectorTags.hpp"
26 #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
28 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
29 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags/Formulation.hpp"
33 #include "Time/Tags.hpp"
35 #include "Utilities/Algorithm.hpp"
40 #include "Utilities/TaggedTuple.hpp"
44 template <
typename... BoundaryCorrectionTags,
typename... Tags,
45 typename BoundaryCorrection>
46 void boundary_correction(
47 const gsl::not_null<Variables<tmpl::list<BoundaryCorrectionTags...>>*>
48 boundary_corrections_on_mortar,
49 const Variables<tmpl::list<Tags...>>& local_boundary_data,
50 const Variables<tmpl::list<Tags...>>& neighbor_boundary_data,
51 const BoundaryCorrection& boundary_correction,
53 boundary_correction.dg_boundary_terms(
55 &get<BoundaryCorrectionTags>(*boundary_corrections_on_mortar))...,
56 get<Tags>(local_boundary_data)..., get<Tags>(neighbor_boundary_data)...,
79 template <
typename Metavariables>
83 Metavariables::volume_dim>>;
84 using const_global_cache_tags = tmpl::list<
88 template <
typename DbTagsList,
typename... InboxTags,
typename ArrayIndex,
89 typename ActionList,
typename ParallelComponent>
93 const ArrayIndex& , ActionList ,
94 const ParallelComponent*
const ) noexcept {
95 constexpr
size_t volume_dim = Metavariables::system::volume_dim;
98 .number_of_neighbors() == 0)) {
100 return {std::move(box)};
103 if constexpr (Metavariables::local_time_stepping) {
110 return {std::move(box)};
113 template <
typename DbTags,
typename... InboxTags,
typename ArrayIndex>
114 static bool is_ready(
const db::DataBox<DbTags>& box,
117 const ArrayIndex& ) noexcept;
120 template <
typename DbTagsList,
typename... InboxTags>
121 static void complete_time_step(
124 template <
typename DbTagsList,
typename... InboxTags>
125 static void receive_global_time_stepping(
129 template <
typename DbTagsList,
typename... InboxTags>
130 static void receive_local_time_stepping(
135 template <
typename Metavariables>
136 template <
typename DbTagsList,
typename... InboxTags>
140 constexpr
size_t volume_dim = Metavariables::system::volume_dim;
150 boost::hash<Key>>>& inbox =
152 volume_dim>>(*inboxes);
153 const auto& temporal_id = get<::Tags::TimeStepId>(*box);
154 const auto received_temporal_id_and_data = inbox.find(temporal_id);
155 db::mutate<evolution::dg::Tags::MortarData<volume_dim>,
158 [&received_temporal_id_and_data](
164 mortar_next_time_step_id) noexcept {
165 for (
auto& received_mortar_data :
166 received_temporal_id_and_data->second) {
167 const auto& mortar_id = received_mortar_data.first;
168 ASSERT(mortar_next_time_step_id->at(mortar_id) ==
169 received_temporal_id_and_data->first,
170 "Expected to receive mortar data on mortar "
171 << mortar_id <<
" at time "
172 << mortar_next_time_step_id->at(mortar_id)
173 <<
" but actually received at time "
174 << received_temporal_id_and_data->first);
175 mortar_next_time_step_id->at(mortar_id) =
176 std::get<3>(received_mortar_data.second);
177 mortar_data->at(mortar_id).insert_neighbor_mortar_data(
178 received_temporal_id_and_data->first,
179 std::get<0>(received_mortar_data.second),
180 std::move(*std::get<2>(received_mortar_data.second)));
183 inbox.erase(received_temporal_id_and_data);
186 template <
typename Metavariables>
187 template <
typename DbTagsList,
typename... InboxTags>
188 void ApplyBoundaryCorrections<Metavariables>::receive_local_time_stepping(
191 constexpr
size_t volume_dim = Metavariables::system::volume_dim;
193 using variables_tag =
typename Metavariables::system::variables_tag;
208 boost::hash<Key>>>& inbox =
210 volume_dim>>(*inboxes);
211 const auto& local_next_temporal_id =
212 db::get<::Tags::Next<::Tags::TimeStepId>>(*box);
215 volume_dim,
typename dt_variables_tag::type>,
217 box, [&inbox, &local_next_temporal_id](
223 typename dt_variables_tag::type>,
225 boundary_data_history,
228 mortar_next_time_step_id) noexcept {
230 for (
auto received_data = inbox.begin();
231 received_data != inbox.end() and
232 received_data->first < local_next_temporal_id;
233 received_data = inbox.erase(received_data)) {
234 const auto& receive_temporal_id = received_data->first;
236 for (
auto& received_mortar_data : received_data->second) {
237 const auto& mortar_id = received_mortar_data.first;
238 MortarData<Metavariables::volume_dim> neighbor_mortar_data{};
243 ASSERT(std::get<2>(received_mortar_data.second).has_value(),
244 "Did not receive boundary correction data from the "
245 "neighbor\nMortarId: "
246 << mortar_id <<
"\nTimeStepId: " << receive_temporal_id);
248 mortar_next_time_step_id->at(mortar_id) == receive_temporal_id,
249 "Expected to receive mortar data on mortar "
250 << mortar_id <<
" at time "
251 << mortar_next_time_step_id->at(mortar_id)
252 <<
" but actually received at time "
253 << receive_temporal_id);
254 mortar_next_time_step_id->at(mortar_id) =
255 std::get<3>(received_mortar_data.second);
256 neighbor_mortar_data.insert_neighbor_mortar_data(
257 receive_temporal_id, std::get<0>(received_mortar_data.second),
258 std::move(*std::get<2>(received_mortar_data.second)));
259 boundary_data_history->at(mortar_id).remote_insert(
260 receive_temporal_id, std::move(neighbor_mortar_data));
266 template <
typename Metavariables>
267 template <
typename DbTagsList,
typename... InboxTags>
268 void ApplyBoundaryCorrections<Metavariables>::complete_time_step(
269 const gsl::not_null<db::DataBox<DbTagsList>*> box) noexcept {
270 constexpr
size_t volume_dim = Metavariables::system::volume_dim;
272 using variables_tag =
typename Metavariables::system::variables_tag;
273 using variables_tags =
typename variables_tag::tags_list;
278 db::get<domain::Tags::Mesh<volume_dim>>(*box);
281 make_array<volume_dim>(volume_mesh.
quadrature(0)),
282 "Must have isotropic quadrature, but got volume mesh: " << volume_mesh);
283 const bool using_gauss_lobatto_points =
284 volume_mesh.
quadrature(0) == Spectral::Quadrature::GaussLobatto;
286 const bool local_time_stepping = Metavariables::local_time_stepping;
289 if (not local_time_stepping and not using_gauss_lobatto_points) {
290 get(volume_det_inv_jacobian)
295 get(volume_det_jacobian) = 1.0 /
get(volume_det_inv_jacobian);
298 const auto& mortar_meshes = db::get<Tags::MortarMesh<volume_dim>>(*box);
299 const auto& mortar_sizes = db::get<Tags::MortarSize<volume_dim>>(*box);
301 db::get<::dg::Tags::Formulation>(*box);
307 face_normal_covector_and_magnitude =
308 db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<volume_dim>>(
311 const TimeDelta& time_step = db::get<::Tags::TimeStep>(*box);
312 const auto& time_stepper =
313 db::get<typename Metavariables::time_stepper_tag>(*box);
315 const auto compute_and_lift_boundary_corrections =
316 [&dg_formulation, &face_normal_covector_and_magnitude,
317 local_time_stepping, &mortar_meshes, &mortar_sizes, &time_step,
318 &time_stepper, using_gauss_lobatto_points, &volume_det_jacobian,
319 &volume_det_inv_jacobian, &volume_mesh](
320 const auto dt_variables_ptr,
const auto variables_ptr,
321 const auto mortar_data_ptr,
const auto mortar_data_history_ptr,
322 const auto& boundary_correction) noexcept {
323 using mortar_tags_list =
typename std::decay_t<decltype(
324 boundary_correction)>::dg_package_field_tags;
326 Variables<db::wrap_tags_in<::Tags::dt, variables_tags>>
327 dt_boundary_correction_on_mortar{};
328 Variables<db::wrap_tags_in<::Tags::dt, variables_tags>>
329 dt_boundary_correction_projected_onto_face{};
335 const auto compute_correction_coupling =
336 [&boundary_correction, dg_formulation,
337 &dt_boundary_correction_on_mortar,
338 &dt_boundary_correction_projected_onto_face, &dt_variables_ptr,
339 &face_normal_covector_and_magnitude, local_time_stepping,
340 &mortar_id_ptr, &mortar_meshes, &mortar_sizes,
341 using_gauss_lobatto_points, &volume_det_jacobian,
342 &volume_det_inv_jacobian, &volume_mesh](
343 const MortarData<volume_dim>& local_mortar_data,
344 const MortarData<volume_dim>& neighbor_mortar_data) noexcept
347 (void)local_time_stepping;
348 ASSERT(mortar_id_ptr !=
nullptr,
349 "Mortar ID pointer should never be nullptr, it must not have "
350 "been set before invoking the lambda.");
351 if (local_time_stepping and not using_gauss_lobatto_points) {
357 local_mortar_data.get_local_volume_det_inv_jacobian(
359 get(volume_det_jacobian) = 1.0 /
get(volume_det_inv_jacobian);
361 const auto& mortar_id = *mortar_id_ptr;
362 const auto& direction = mortar_id.first;
363 const auto&
mortar_mesh = mortar_meshes.at(mortar_id);
368 local_mesh_and_data = *local_mortar_data.local_mortar_data();
370 neighbor_mesh_and_data =
371 *neighbor_mortar_data.neighbor_mortar_data();
372 Variables<mortar_tags_list> local_data_on_mortar{
374 Variables<mortar_tags_list> neighbor_data_on_mortar{
376 std::copy(std::get<1>(local_mesh_and_data).
begin(),
377 std::get<1>(local_mesh_and_data).
end(),
378 local_data_on_mortar.data());
379 std::copy(std::get<1>(neighbor_mesh_and_data).
begin(),
380 std::get<1>(neighbor_mesh_and_data).
end(),
381 neighbor_data_on_mortar.data());
390 if (dt_boundary_correction_on_mortar.number_of_grid_points() !=
392 dt_boundary_correction_on_mortar.initialize(
396 detail::boundary_correction(
398 local_data_on_mortar, neighbor_data_on_mortar,
399 boundary_correction, dg_formulation);
402 mortar_sizes.at(mortar_id);
403 const Mesh<volume_dim - 1> face_mesh =
404 volume_mesh.
slice_away(direction.dimension());
406 auto& dt_boundary_correction =
407 [&dt_boundary_correction_on_mortar,
408 &dt_boundary_correction_projected_onto_face, &face_mesh,
412 dt_boundary_correction_projected_onto_face =
416 return dt_boundary_correction_projected_onto_face;
418 return dt_boundary_correction_on_mortar;
422 if (local_time_stepping) {
423 local_mortar_data.get_local_face_normal_magnitude(
424 &magnitude_of_face_normal);
426 ASSERT(face_normal_covector_and_magnitude.count(direction) == 1 and
427 face_normal_covector_and_magnitude.at(direction)
429 "Face normal covector and magnitude not set in "
432 get(magnitude_of_face_normal)
434 get(get<evolution::dg::Tags::MagnitudeOfNormal>(
435 *face_normal_covector_and_magnitude.at(direction))))));
438 if (using_gauss_lobatto_points) {
442 volume_mesh.
extents(direction.dimension()),
443 magnitude_of_face_normal);
444 if (local_time_stepping) {
445 return dt_boundary_correction;
449 dt_variables_ptr, dt_boundary_correction,
450 volume_mesh.
extents(), direction.dimension(),
466 if (local_time_stepping) {
468 "For local time stepping the volume determinant of the "
469 "inverse Jacobian has not been set.");
472 local_mortar_data.get_local_face_det_jacobian(
475 Variables<db::wrap_tags_in<::Tags::dt, variables_tags>>
476 volume_dt_correction{
477 dt_variables_ptr->number_of_grid_points(), 0.0};
479 make_not_null(&volume_dt_correction), volume_det_inv_jacobian,
480 volume_mesh, direction, dt_boundary_correction,
481 magnitude_of_face_normal, face_det_jacobian);
482 return volume_dt_correction;
487 face_mesh.number_of_grid_points()};
489 auto interpolation_matrices =
490 make_array<volume_dim>(std::cref(
identity));
494 gsl::at(interpolation_matrices, direction.dimension()) =
495 direction.side() == Side::Upper ? matrices.second
498 interpolation_matrices,
get(volume_det_jacobian),
502 dt_variables_ptr, volume_det_inv_jacobian, volume_mesh,
503 direction, dt_boundary_correction, magnitude_of_face_normal,
509 not local_time_stepping,
510 "Must return lifted data when using local time stepping. It's "
511 "likely a missing return statement or incorrect if-else logical "
512 "related to whether or not we are doing local time stepping.");
516 if constexpr (Metavariables::local_time_stepping) {
517 (void)mortar_data_ptr;
519 for (
auto& mortar_id_and_data : *mortar_data_history_ptr) {
520 const auto& mortar_id = mortar_id_and_data.first;
521 auto& mortar_data_history = mortar_id_and_data.second;
522 const auto& direction = mortar_id.first;
526 "Cannot impose boundary conditions on external boundary in "
529 <<
" in the ApplyBoundaryCorrections action. Boundary "
530 "conditions are applied in the ComputeTimeDerivative "
532 "instead. You may have unintentionally added external "
533 "mortars in one of the initialization actions.");
535 mortar_id_ptr = &mortar_id;
536 auto lifted_volume_data = time_stepper.compute_boundary_delta(
537 compute_correction_coupling,
540 if (using_gauss_lobatto_points) {
543 variables_ptr, lifted_volume_data, volume_mesh.
extents(),
544 direction.dimension(),
547 *variables_ptr += lifted_volume_data;
554 (void)mortar_data_history_ptr;
556 for (
auto& mortar_id_and_data : *mortar_data_ptr) {
559 const auto& mortar_id = mortar_id_and_data.first;
560 const auto& mortar_data = mortar_id_and_data.second;
561 const auto& direction = mortar_id.first;
565 "Cannot impose boundary conditions on external boundary in "
568 <<
" in the ApplyBoundaryCorrections action. Boundary "
569 "conditions are applied in the ComputeTimeDerivative "
571 "instead. You may have unintentionally added external "
572 "mortars in one of the initialization actions.");
574 mortar_id_ptr = &mortar_id;
575 compute_correction_coupling(mortar_data, mortar_data);
582 const auto& boundary_correction =
db::get<
585 using derived_boundary_corrections =
586 typename std::decay_t<decltype(boundary_correction)>::creatable_classes;
590 "All createable classes for boundary corrections must be marked "
592 tmpl::for_each<derived_boundary_corrections>(
593 [&boundary_correction, &box, &compute_and_lift_boundary_corrections](
594 auto derived_correction_v) noexcept {
595 using DerivedCorrection =
596 tmpl::type_from<decltype(derived_correction_v)>;
597 if (
typeid(boundary_correction) ==
typeid(DerivedCorrection)) {
604 volume_dim,
typename dt_variables_tag::type>>(
605 box, compute_and_lift_boundary_corrections,
606 dynamic_cast<const DerivedCorrection&
>(boundary_correction));
611 template <
typename Metavariables>
612 template <
typename DbTags,
typename... InboxTags,
typename ArrayIndex>
613 bool ApplyBoundaryCorrections<Metavariables>::is_ready(
614 const db::DataBox<DbTags>& box,
617 const ArrayIndex& ) noexcept {
618 static constexpr
size_t volume_dim = Metavariables::volume_dim;
620 db::get<domain::Tags::Element<volume_dim>>(box);
626 if constexpr (not Metavariables::local_time_stepping) {
627 const TimeStepId& temporal_id = get<::Tags::TimeStepId>(box);
630 Metavariables::volume_dim>>(inboxes);
631 const auto received_temporal_id_and_data = inbox.find(temporal_id);
632 if (received_temporal_id_and_data == inbox.end()) {
635 const auto& received_neighbor_data = received_temporal_id_and_data->second;
636 for (
const auto& [direction, neighbors] : element.
neighbors()) {
637 for (
const auto& neighbor : neighbors) {
638 const auto neighbor_received =
639 received_neighbor_data.find(
std::pair{direction, neighbor});
640 if (neighbor_received == received_neighbor_data.end()) {
648 Metavariables::volume_dim>>(inboxes);
650 const auto& local_next_temporal_id =
651 db::get<::Tags::Next<::Tags::TimeStepId>>(box);
652 const auto& mortars_next_temporal_id =
db::get<
656 for (
const auto& [mortar_id, mortar_next_temporal_id] :
657 mortars_next_temporal_id) {
663 auto next_temporal_id = mortar_next_temporal_id;
664 while (next_temporal_id < local_next_temporal_id) {
665 const auto temporal_received = inbox.find(next_temporal_id);
666 if (temporal_received == inbox.end()) {
669 const auto mortar_received = temporal_received->second.find(mortar_id);
670 if (mortar_received == temporal_received->second.end()) {
673 next_temporal_id = std::get<3>(mortar_received->second);