15 #include "DataStructures/DataBox/PrefixHelpers.hpp"
17 #include "DataStructures/DataVector.hpp"
18 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
22 #include "Domain/BoundaryConditions/None.hpp"
23 #include "Domain/BoundaryConditions/Periodic.hpp"
26 #include "Domain/ElementMap.hpp"
27 #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
28 #include "Domain/FunctionsOfTime/Tags.hpp"
33 #include "Domain/TagsTimeDependent.hpp"
34 #include "Evolution/BoundaryConditions/Type.hpp"
35 #include "Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivativeHelpers.hpp"
36 #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
37 #include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp"
38 #include "Evolution/DiscontinuousGalerkin/InterpolateFromBoundary.hpp"
39 #include "Evolution/DiscontinuousGalerkin/LiftFromBoundary.hpp"
40 #include "Evolution/DiscontinuousGalerkin/ProjectToBoundary.hpp"
41 #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
43 #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags/Formulation.hpp"
47 #include "Time/Tags.hpp"
53 namespace evolution::dg::Actions::detail {
54 template <
typename BoundaryConditionHelper,
typename AllTagsOnFaceList,
55 typename... TagsFromFace,
typename... VolumeArgs>
57 BoundaryConditionHelper& boundary_condition_helper,
58 const Variables<AllTagsOnFaceList>& fields_on_interior_face,
59 tmpl::list<TagsFromFace...> ,
60 const VolumeArgs&... volume_args) noexcept {
61 return boundary_condition_helper(
62 get<TagsFromFace>(fields_on_interior_face)..., volume_args...);
65 template <
typename System,
size_t Dim,
typename DbTagsList,
66 typename BoundaryCorrection,
typename BoundaryCondition,
67 typename... EvolvedVariablesTags,
typename... PackageDataVolumeTags,
68 typename... BoundaryConditionVolumeTags,
typename... PackageFieldTags,
69 typename... BoundaryCorrectionPackagedDataInputTags>
70 void apply_boundary_condition_on_face(
72 [[maybe_unused]]
const BoundaryCorrection& boundary_correction,
73 const BoundaryCondition& boundary_condition,
75 [[maybe_unused]]
const Variables<tmpl::list<EvolvedVariablesTags...>>&
77 [[maybe_unused]]
const Variables<
80 [[maybe_unused]]
const Variables<
83 [[maybe_unused]]
const Variables<
84 typename System::compute_volume_time_derivative_terms::temporary_tags>&
86 [[maybe_unused]]
const Variables<
87 detail::get_primitive_vars_tags_from_system<System>>*
const
88 volume_primitive_variables,
91 [[maybe_unused]] const ::ElementMap<Dim, Frame::Grid>& logical_to_grid_map,
94 [[maybe_unused]]
const double time,
98 const std::optional<tnsr::I<DataVector, Dim>>& volume_mesh_velocity,
99 const InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>&
100 volume_inverse_jacobian,
102 tmpl::list<PackageDataVolumeTags...> ,
103 tmpl::list<PackageFieldTags...> ,
104 tmpl::list<BoundaryCorrectionPackagedDataInputTags...> ,
105 tmpl::list<BoundaryConditionVolumeTags...> ) noexcept {
106 using variables_tag =
typename System::variables_tag;
107 using variables_tags =
typename variables_tag::tags_list;
108 using flux_variables =
typename System::flux_variables;
112 const Mesh<Dim - 1> face_mesh = volume_mesh.
slice_away(direction.dimension());
113 const size_t number_of_points_on_face = face_mesh.number_of_grid_points();
126 constexpr
bool uses_ghost_condition =
127 BoundaryCondition::bc_type ==
128 evolution::BoundaryConditions::Type::Ghost or
129 BoundaryCondition::bc_type ==
130 evolution::BoundaryConditions::Type::GhostAndTimeDerivative;
131 constexpr
bool uses_time_derivative_condition =
132 BoundaryCondition::bc_type ==
133 evolution::BoundaryConditions::Type::TimeDerivative or
134 BoundaryCondition::bc_type ==
135 evolution::BoundaryConditions::Type::GhostAndTimeDerivative;
136 constexpr
bool needs_coordinates = tmpl::list_contains_v<
137 typename BoundaryCondition::dg_interior_temporary_tags,
141 using inverse_spatial_metric_list =
142 detail::inverse_spatial_metric_tag<System>;
143 constexpr
bool has_inv_spatial_metric =
144 detail::has_inverse_spatial_metric_tag_v<System>;
147 using bcondition_interior_temp_tags =
148 typename BoundaryCondition::dg_interior_temporary_tags;
149 using bcondition_interior_prim_tags =
150 detail::boundary_condition_primitive_tags<
151 System::has_primitive_and_conservative_vars, BoundaryCondition>;
152 using bcondition_interior_evolved_vars_tags =
153 typename BoundaryCondition::dg_interior_evolved_variables_tags;
154 using bcondition_interior_dt_evolved_vars_tags =
155 detail::get_dt_vars_from_boundary_condition<BoundaryCondition>;
156 using bcondition_interior_deriv_evolved_vars_tags =
157 detail::get_deriv_vars_from_boundary_condition<BoundaryCondition>;
158 using bcondition_interior_tags = tmpl::append<
159 tmpl::conditional_t<has_inv_spatial_metric,
160 tmpl::list<detail::NormalVector<Dim>>, tmpl::list<>>,
161 bcondition_interior_evolved_vars_tags, bcondition_interior_prim_tags,
162 bcondition_interior_temp_tags, bcondition_interior_dt_evolved_vars_tags,
163 bcondition_interior_deriv_evolved_vars_tags>;
166 using correction_temp_tags = tmpl::conditional_t<
167 uses_ghost_condition,
168 typename BoundaryCorrection::dg_package_data_temporary_tags,
170 using correction_prim_tags = tmpl::conditional_t<
171 uses_ghost_condition,
172 detail::boundary_correction_primitive_tags<
173 System::has_primitive_and_conservative_vars, BoundaryCorrection>,
175 using correction_evolved_vars_tags =
176 tmpl::conditional_t<uses_ghost_condition,
177 typename System::variables_tag::tags_list,
184 using interior_temp_tags = tmpl::remove_duplicates<
185 tmpl::append<bcondition_interior_temp_tags, correction_temp_tags>>;
186 using interior_prim_tags = tmpl::remove_duplicates<
187 tmpl::append<bcondition_interior_prim_tags, correction_prim_tags>>;
188 using interior_evolved_vars_tags = tmpl::remove_duplicates<tmpl::append<
189 correction_evolved_vars_tags, bcondition_interior_evolved_vars_tags>>;
195 tmpl::conditional_t<uses_ghost_condition,
199 using tags_on_interior_face = tmpl::append<
200 fluxes_tags, interior_temp_tags, interior_prim_tags,
201 interior_evolved_vars_tags, bcondition_interior_dt_evolved_vars_tags,
202 bcondition_interior_deriv_evolved_vars_tags, inverse_spatial_metric_list,
203 tmpl::list<detail::OneOverNormalVectorMagnitude,
204 detail::NormalVector<Dim>>>;
206 Variables<tags_on_interior_face> interior_face_fields{
207 number_of_points_on_face};
220 if constexpr (uses_ghost_condition) {
222 volume_evolved_vars, volume_mesh,
225 project_tensors_to_boundary<interior_evolved_vars_tags>(
226 make_not_null(&interior_face_fields), volume_evolved_vars, volume_mesh,
229 if constexpr (tmpl::size<fluxes_tags>::value != 0) {
231 volume_fluxes, volume_mesh, direction);
235 using temp_tags_no_coordinates =
236 tmpl::remove<interior_temp_tags,
238 if constexpr (tmpl::size<tmpl::append<
239 temp_tags_no_coordinates,
240 detail::inverse_spatial_metric_tag<System>>>::value != 0) {
242 temp_tags_no_coordinates, detail::inverse_spatial_metric_tag<System>>>(
243 make_not_null(&interior_face_fields), volume_temporaries, volume_mesh,
246 if constexpr (System::has_primitive_and_conservative_vars and
247 tmpl::size<interior_prim_tags>::value != 0) {
248 ASSERT(volume_primitive_variables !=
nullptr,
249 "The volume primitive variables are not set even though the "
250 "system has primitive variables.");
251 project_tensors_to_boundary<interior_prim_tags>(
252 make_not_null(&interior_face_fields), *volume_primitive_variables,
253 volume_mesh, direction);
255 (void)volume_primitive_variables;
257 if constexpr (tmpl::size<
258 bcondition_interior_deriv_evolved_vars_tags>::value != 0) {
259 project_tensors_to_boundary<bcondition_interior_deriv_evolved_vars_tags>(
260 make_not_null(&interior_face_fields), partial_derivs, volume_mesh,
263 if constexpr (tmpl::size<bcondition_interior_dt_evolved_vars_tags>::value !=
265 project_tensors_to_boundary<bcondition_interior_dt_evolved_vars_tags>(
266 make_not_null(&interior_face_fields), db::get<dt_variables_tag>(*box),
267 volume_mesh, direction);
271 if (volume_mesh_velocity.has_value()) {
272 face_mesh_velocity = tnsr::I<DataVector, Dim>{number_of_points_on_face};
274 *volume_mesh_velocity, volume_mesh, direction);
279 const auto normalize_normal_vectors =
280 [&direction, mesh_is_moving = not moving_mesh_map.is_identity(),
281 number_of_points_on_face, &volume_inverse_jacobian,
282 &volume_mesh](
const auto normal_covector_magnitude_in_direction_ptr,
283 auto fields_on_face_ptr) noexcept {
284 if (
auto& normal_covector_quantity =
285 *normal_covector_magnitude_in_direction_ptr;
286 has_inv_spatial_metric or mesh_is_moving or
287 not normal_covector_quantity.has_value()) {
288 if (not normal_covector_quantity.has_value()) {
289 normal_covector_quantity =
292 number_of_points_on_face};
294 tnsr::i<DataVector, Dim> volume_unnormalized_normal_covector{};
296 for (
size_t inertial_index = 0; inertial_index < Dim;
298 volume_unnormalized_normal_covector.get(inertial_index)
301 volume_inverse_jacobian
302 .
get(direction.dimension(), inertial_index)
304 volume_mesh.number_of_grid_points());
308 *normal_covector_quantity)),
309 volume_unnormalized_normal_covector, volume_mesh, direction);
311 if (
const double sign = direction.sign(); sign != 1.0) {
312 for (
auto& normal_covector_component :
314 *normal_covector_quantity)) {
315 normal_covector_component *= sign;
319 detail::unit_normal_vector_and_covector_and_magnitude_impl<System>(
321 *normal_covector_quantity)),
323 *normal_covector_quantity)),
326 *normal_covector_quantity));
330 db::mutate<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>>(
331 box, [&direction, &interior_face_fields, &normalize_normal_vectors](
332 const auto normal_covector_and_magnitude_ptr) noexcept {
333 normalize_normal_vectors(
334 make_not_null(&normal_covector_and_magnitude_ptr->at(direction)),
338 const tnsr::i<DataVector, Dim, Frame::Inertial>& interior_normal_covector =
339 get<evolution::dg::Tags::NormalCovector<Dim>>(
340 *db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>>(*box)
343 if constexpr (needs_coordinates) {
345 get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(interior_face_fields) =
347 face_mesh, direction)),
348 time, functions_of_time);
351 if constexpr (BoundaryCondition::bc_type ==
352 evolution::BoundaryConditions::Type::Outflow) {
356 const auto apply_bc = [&boundary_condition, &face_mesh_velocity,
357 &interior_normal_covector](
358 const auto&... face_and_volume_args) noexcept {
359 return boundary_condition.dg_outflow(face_mesh_velocity,
360 interior_normal_covector,
361 face_and_volume_args...);
364 apply_boundary_condition_impl(
365 apply_bc, interior_face_fields, bcondition_interior_tags{},
366 db::get<BoundaryConditionVolumeTags>(*box)...);
367 if (error_message.has_value()) {
368 ERROR(*error_message <<
"\n\nIn element:" << element.id()
369 <<
"\nIn direction: " << direction);
378 Variables<dt_variables_tags> dt_time_derivative_correction{};
379 if constexpr (uses_time_derivative_condition) {
380 dt_time_derivative_correction.initialize(number_of_points_on_face);
381 auto apply_bc = [&boundary_condition, &dt_time_derivative_correction,
382 &face_mesh_velocity, &interior_normal_covector](
383 const auto&... interior_face_and_volume_args) {
384 return boundary_condition.dg_time_derivative(
386 dt_time_derivative_correction))...,
387 face_mesh_velocity, interior_normal_covector,
388 interior_face_and_volume_args...);
391 apply_boundary_condition_impl(
392 apply_bc, interior_face_fields, bcondition_interior_tags{},
393 db::get<BoundaryConditionVolumeTags>(*box)...);
394 if (error_message.has_value()) {
395 ERROR(*error_message <<
"\n\nIn element:" << element.id()
396 <<
"\nIn direction: " << direction);
399 (void)dt_time_derivative_correction;
404 using tags_on_exterior_face =
405 tmpl::append<variables_tags, fluxes_tags, correction_temp_tags,
406 correction_prim_tags, inverse_spatial_metric_list,
407 tmpl::list<detail::OneOverNormalVectorMagnitude,
408 detail::NormalVector<Dim>,
410 Variables<tags_on_exterior_face> exterior_face_fields{
411 number_of_points_on_face};
413 if constexpr (uses_ghost_condition) {
414 using mortar_tags_list = tmpl::list<PackageFieldTags...>;
415 using dg_package_data_projected_tags =
416 tmpl::append<variables_tags, fluxes_tags, correction_temp_tags,
417 correction_prim_tags>;
419 Variables<mortar_tags_list> internal_packaged_data{
420 number_of_points_on_face};
421 const double max_abs_char_speed_on_face = detail::dg_package_data<System>(
423 interior_face_fields, interior_normal_covector, face_mesh_velocity,
424 dg_package_data_projected_tags{},
425 db::get<PackageDataVolumeTags>(*box)...);
426 (void)max_abs_char_speed_on_face;
445 auto apply_bc = [&boundary_condition, &exterior_face_fields,
446 &face_mesh_velocity, &interior_normal_covector](
447 const auto&... interior_face_and_volume_args) {
448 if constexpr (has_inv_spatial_metric) {
449 return boundary_condition.dg_ghost(
451 exterior_face_fields))...,
453 &
get<tmpl::front<detail::inverse_spatial_metric_tag<System>>>(
454 exterior_face_fields)),
455 face_mesh_velocity, interior_normal_covector,
456 interior_face_and_volume_args...);
458 return boundary_condition.dg_ghost(
460 exterior_face_fields))...,
461 face_mesh_velocity, interior_normal_covector,
462 interior_face_and_volume_args...);
466 apply_boundary_condition_impl(
467 apply_bc, interior_face_fields, bcondition_interior_tags{},
468 db::get<BoundaryConditionVolumeTags>(*box)...);
469 if (error_message.has_value()) {
470 ERROR(*error_message <<
"\n\nIn element:" << element.id()
471 <<
"\nIn direction: " << direction);
474 if (face_mesh_velocity.has_value()) {
475 tmpl::for_each<flux_variables>(
476 [&face_mesh_velocity, &exterior_face_fields](
auto tag_v) noexcept {
478 using var_tag =
typename decltype(tag_v)::type;
482 auto& flux_var = get<flux_var_tag>(exterior_face_fields);
483 const auto& var = get<var_tag>(exterior_face_fields);
484 const auto& mesh_velocity = *face_mesh_velocity;
486 for (
size_t flux_var_storage_index = 0;
487 flux_var_storage_index < flux_var.size();
488 ++flux_var_storage_index) {
491 const auto flux_var_tensor_index =
492 flux_var.get_tensor_index(flux_var_storage_index);
495 const auto var_tensor_index =
498 const size_t flux_index =
gsl::at(flux_var_tensor_index, 0);
501 flux_var[flux_var_storage_index] -=
502 var.get(var_tensor_index) * mesh_velocity.get(flux_index);
510 for (
size_t i = 0; i < Dim; ++i) {
511 get<evolution::dg::Tags::NormalCovector<Dim>>(exterior_face_fields)
512 .
get(i) = -interior_normal_covector.get(i);
514 if constexpr (has_inv_spatial_metric) {
515 const tnsr::II<DataVector, Dim, Frame::Inertial>& inv_spatial_metric =
516 get<tmpl::front<inverse_spatial_metric_list>>(exterior_face_fields);
517 tnsr::i<DataVector, Dim, Frame::Inertial>& exterior_normal_covector =
518 get<evolution::dg::Tags::NormalCovector<Dim>>(exterior_face_fields);
519 tnsr::I<DataVector, Dim, Frame::Inertial>& exterior_normal_vector =
520 get<detail::NormalVector<Dim>>(exterior_face_fields);
525 for (
size_t i = 0; i < Dim; ++i) {
526 exterior_normal_vector.get(i) =
527 get<0>(exterior_normal_covector) * inv_spatial_metric.get(i, 0);
528 for (
size_t j = 1; j < Dim; ++j) {
529 exterior_normal_vector.get(i) +=
530 exterior_normal_covector.get(j) * inv_spatial_metric.get(i, j);
537 get<detail::OneOverNormalVectorMagnitude>(exterior_face_fields);
539 exterior_normal_vector);
541 for (
size_t i = 0; i < Dim; ++i) {
548 Variables<mortar_tags_list> external_packaged_data{
549 number_of_points_on_face};
550 detail::dg_package_data<System>(
552 exterior_face_fields,
554 face_mesh_velocity, dg_package_data_projected_tags{},
555 db::get<PackageDataVolumeTags>(*box)...);
557 Variables<dt_variables_tags> boundary_corrections_on_face{
558 number_of_points_on_face};
561 boundary_correction.dg_boundary_terms(
563 boundary_corrections_on_face))...,
564 get<PackageFieldTags>(internal_packaged_data)...,
565 get<PackageFieldTags>(external_packaged_data)..., dg_formulation);
568 const auto& magnitude_of_interior_face_normal =
569 get<evolution::dg::Tags::MagnitudeOfNormal>(
572 if (volume_mesh.quadrature(0) == Spectral::Quadrature::GaussLobatto) {
576 volume_mesh.extents(direction.dimension()),
577 magnitude_of_interior_face_normal);
580 db::mutate<dt_variables_tag>(
581 box, [&direction, &boundary_corrections_on_face,
582 &volume_mesh](
const auto dt_variables_ptr) noexcept {
584 dt_variables_ptr, boundary_corrections_on_face,
585 volume_mesh.extents(), direction.dimension(),
594 const DataVector volume_det_jacobian = 1.0 /
get(volume_det_inv_jacobian);
600 auto interpolation_matrices = make_array<Dim>(std::cref(
identity));
603 volume_mesh.slice_through(direction.dimension()));
604 gsl::at(interpolation_matrices, direction.dimension()) =
605 direction.side() == Side::Upper ? matrices.second : matrices.first;
607 interpolation_matrices, volume_det_jacobian,
608 volume_mesh.extents());
610 db::mutate<dt_variables_tag>(
611 box, [&direction, &boundary_corrections_on_face, &face_det_jacobian,
612 &magnitude_of_interior_face_normal, &volume_det_inv_jacobian,
613 &volume_mesh](
const auto dt_variables_ptr) noexcept {
615 dt_variables_ptr, volume_det_inv_jacobian, volume_mesh,
616 direction, boundary_corrections_on_face,
617 magnitude_of_interior_face_normal, face_det_jacobian);
622 if constexpr (uses_time_derivative_condition) {
623 if (volume_mesh.quadrature(0) == Spectral::Quadrature::GaussLobatto) {
624 db::mutate<dt_variables_tag>(
625 box, [&direction, &dt_time_derivative_correction,
626 &volume_mesh](
const auto dt_variables_ptr) noexcept {
628 dt_variables_ptr, dt_time_derivative_correction,
629 volume_mesh.extents(), direction.dimension(),
633 db::mutate<dt_variables_tag>(
634 box, [&direction, &dt_time_derivative_correction,
635 &volume_mesh](
const auto dt_variables_ptr) noexcept {
638 dt_time_derivative_correction);
653 template <
typename System,
size_t Dim,
typename DbTagsList,
654 typename BoundaryCorrection>
655 void apply_boundary_conditions_on_all_external_faces(
657 const BoundaryCorrection& boundary_correction,
659 typename System::compute_volume_time_derivative_terms::temporary_tags>&
667 const Variables<detail::get_primitive_vars_tags_from_system<System>>*
const
668 primitive_vars) noexcept {
669 using derived_boundary_conditions = tmpl::remove_if<
670 typename System::boundary_conditions_base::creatable_classes,
676 using variables_tag =
typename System::variables_tag;
677 using flux_variables =
typename System::flux_variables;
681 const Element<Dim>& element = db::get<domain::Tags::Element<Dim>>(*box);
684 if (number_of_boundaries_left == 0) {
688 const auto& external_boundary_conditions =
689 db::get<domain::Tags::Domain<Dim>>(*box)
690 .blocks()[element.
id().block_id()]
691 .external_boundary_conditions();
692 tmpl::for_each<derived_boundary_conditions>(
693 [&boundary_correction, &box, &element, &external_boundary_conditions,
694 &number_of_boundaries_left, &partial_derivs, &primitive_vars,
696 &volume_fluxes](
auto derived_boundary_condition_v) noexcept {
697 using DerivedBoundaryCondition =
698 tmpl::type_from<decltype(derived_boundary_condition_v)>;
700 if (number_of_boundaries_left == 0) {
705 const auto& boundary_condition =
706 *external_boundary_conditions.at(direction);
707 if (
typeid(boundary_condition) ==
typeid(DerivedBoundaryCondition)) {
708 detail::apply_boundary_condition_on_face<System>(
709 box, boundary_correction,
710 dynamic_cast<const DerivedBoundaryCondition&
>(
712 direction, db::get<variables_tag>(*box), volume_fluxes,
713 partial_derivs, temporaries, primitive_vars,
714 db::get<::dg::Tags::Formulation>(*box),
720 db::get<::Tags::Time>(*box),
721 db::get<::domain::Tags::FunctionsOfTime>(*box),
727 typename BoundaryCorrection::dg_package_data_volume_tags{},
728 typename BoundaryCorrection::dg_package_field_tags{},
730 typename variables_tag::tags_list, fluxes_tags,
731 typename BoundaryCorrection::dg_package_data_temporary_tags,
732 typename detail::get_primitive_vars<
733 System::has_primitive_and_conservative_vars>::
734 template f<BoundaryCorrection>>{},
735 typename DerivedBoundaryCondition::dg_gridless_tags{});
736 --number_of_boundaries_left;
738 if (number_of_boundaries_left == 0) {