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 <functional>
8 : #include <memory>
9 : #include <optional>
10 : #include <string>
11 : #include <tuple>
12 : #include <unordered_map>
13 : #include <utility>
14 :
15 : #include "DataStructures/DataBox/DataBox.hpp"
16 : #include "DataStructures/DataBox/MetavariablesTag.hpp"
17 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
18 : #include "DataStructures/DataBox/Prefixes.hpp"
19 : #include "DataStructures/DataVector.hpp"
20 : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
21 : #include "DataStructures/Tensor/Tensor.hpp"
22 : #include "DataStructures/Variables.hpp"
23 : #include "DataStructures/VariablesTag.hpp"
24 : #include "Domain/Block.hpp"
25 : #include "Domain/BoundaryConditions/Cartoon.hpp"
26 : #include "Domain/BoundaryConditions/None.hpp"
27 : #include "Domain/BoundaryConditions/Periodic.hpp"
28 : #include "Domain/CoordinateMaps/CoordinateMap.hpp"
29 : #include "Domain/Creators/Tags/ExternalBoundaryConditions.hpp"
30 : #include "Domain/Domain.hpp"
31 : #include "Domain/ElementMap.hpp"
32 : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
33 : #include "Domain/FunctionsOfTime/Tags.hpp"
34 : #include "Domain/InterfaceLogicalCoordinates.hpp"
35 : #include "Domain/Structure/Direction.hpp"
36 : #include "Domain/Structure/Element.hpp"
37 : #include "Domain/Structure/ElementId.hpp"
38 : #include "Domain/Tags.hpp"
39 : #include "Domain/TagsTimeDependent.hpp"
40 : #include "Evolution/BoundaryConditions/Type.hpp"
41 : #include "Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivativeHelpers.hpp"
42 : #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
43 : #include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp"
44 : #include "Evolution/DiscontinuousGalerkin/BoundaryEvolvedVariables.hpp"
45 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
46 : #include "NumericalAlgorithms/DiscontinuousGalerkin/InterpolateFromBoundary.hpp"
47 : #include "NumericalAlgorithms/DiscontinuousGalerkin/LiftFlux.hpp"
48 : #include "NumericalAlgorithms/DiscontinuousGalerkin/LiftFromBoundary.hpp"
49 : #include "NumericalAlgorithms/DiscontinuousGalerkin/ProjectToBoundary.hpp"
50 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags/Formulation.hpp"
51 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
52 : #include "NumericalAlgorithms/Spectral/BoundaryInterpolationMatrices.hpp"
53 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
54 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
55 : #include "Utilities/ErrorHandling/Assert.hpp"
56 : #include "Utilities/ErrorHandling/Error.hpp"
57 : #include "Utilities/Gsl.hpp"
58 : #include "Utilities/TMPL.hpp"
59 :
60 : /// \cond
61 : namespace Tags {
62 : struct Time;
63 : } // namespace Tags
64 : /// \endcond
65 :
66 : namespace evolution::dg::Actions::detail {
67 : template <typename BoundaryConditionHelper, typename AllTagsOnFaceList,
68 : typename... TagsFromFace, typename... VolumeArgs>
69 : std::optional<std::string> apply_boundary_condition_impl(
70 : BoundaryConditionHelper& boundary_condition_helper,
71 : const Variables<AllTagsOnFaceList>& fields_on_interior_face,
72 : tmpl::list<TagsFromFace...> /*meta*/, const VolumeArgs&... volume_args) {
73 : return boundary_condition_helper(
74 : get<TagsFromFace>(fields_on_interior_face)..., volume_args...);
75 : }
76 :
77 : template <typename System, size_t Dim, bool ComputeAuxiliary = false,
78 : typename DbTagsList, typename BoundaryCorrection,
79 : typename BoundaryCondition, typename... EvolvedVariablesTags,
80 : typename... PackageDataVolumeTags,
81 : typename... BoundaryConditionVolumeTags, typename... PackageFieldTags,
82 : typename... BoundaryTermsVolumeTags,
83 : typename... BoundaryCorrectionPackagedDataInputTags,
84 : typename... AuxiliaryCorrectionTags>
85 : void apply_boundary_condition_on_face(
86 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
87 : [[maybe_unused]] const BoundaryCorrection& boundary_correction,
88 : const BoundaryCondition& boundary_condition,
89 : const Direction<Dim>& direction,
90 : [[maybe_unused]] const Variables<tmpl::list<EvolvedVariablesTags...>>&
91 : volume_evolved_vars,
92 : // The auxiliary variables are projected to the face in both passes
93 : // because dg_ghost is called in the auxiliary pass and may declare
94 : // auxiliary variables in dg_interior_evolved_variables_tags.
95 : [[maybe_unused]] const Variables<
96 : get_auxiliary_variables_or_default_t<System, tmpl::list<>>>* const
97 : volume_auxiliary_variables,
98 : [[maybe_unused]] const Variables<
99 : db::wrap_tags_in<::Tags::Flux, typename System::flux_variables,
100 : tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes,
101 : [[maybe_unused]] const Variables<
102 : db::wrap_tags_in<::Tags::deriv, typename System::gradient_variables,
103 : tmpl::size_t<Dim>, Frame::Inertial>>& partial_derivs,
104 : [[maybe_unused]] const Variables<
105 : typename System::compute_volume_time_derivative_terms::temporary_tags>&
106 : volume_temporaries,
107 : [[maybe_unused]] const Variables<
108 : detail::get_primitive_vars_tags_from_system<System>>* const
109 : volume_primitive_variables,
110 : [[maybe_unused]] const ::dg::Formulation dg_formulation,
111 : const Mesh<Dim>& volume_mesh, [[maybe_unused]] const Element<Dim>& element,
112 : [[maybe_unused]] const ::ElementMap<Dim, Frame::Grid>& logical_to_grid_map,
113 : const domain::CoordinateMapBase<Frame::Grid, Frame::Inertial, Dim>&
114 : moving_mesh_map,
115 : [[maybe_unused]] const double time,
116 : [[maybe_unused]] const std::unordered_map<
117 : std::string, std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>&
118 : functions_of_time,
119 : const std::optional<tnsr::I<DataVector, Dim>>& volume_mesh_velocity,
120 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
121 : Frame::Inertial>& volume_inverse_jacobian,
122 : [[maybe_unused]] const Scalar<DataVector>& volume_det_inv_jacobian,
123 : tmpl::list<PackageDataVolumeTags...> /*meta*/,
124 : tmpl::list<PackageFieldTags...> /*meta*/,
125 : tmpl::list<BoundaryTermsVolumeTags...> /*meta*/,
126 : tmpl::list<BoundaryCorrectionPackagedDataInputTags...> /*meta*/,
127 : tmpl::list<BoundaryConditionVolumeTags...> /*meta*/,
128 : tmpl::list<AuxiliaryCorrectionTags...> /*meta*/) {
129 : using variables_tags = tmpl::list<EvolvedVariablesTags...>;
130 : using variables_tag = ::Tags::Variables<variables_tags>;
131 : using flux_variables = typename System::flux_variables;
132 : using dt_variables_tags = db::wrap_tags_in<::Tags::dt, variables_tags>;
133 : using dt_variables_tag = db::add_tag_prefix<::Tags::dt, variables_tag>;
134 : using auxiliary_variables =
135 : get_auxiliary_variables_or_default_t<System, tmpl::list<>>;
136 : using auxiliary_variables_tag = ::Tags::Variables<auxiliary_variables>;
137 : using tag_to_update =
138 : tmpl::conditional_t<ComputeAuxiliary, auxiliary_variables_tag,
139 : dt_variables_tag>;
140 : using projected_auxiliary_vars_tags =
141 : tmpl::conditional_t<ComputeAuxiliary, tmpl::list<>, auxiliary_variables>;
142 :
143 : constexpr bool bc_evolves_boundary_variables =
144 : evolution::dg::evolves_boundary_variables_v<BoundaryCondition>;
145 : static_assert(
146 : not bc_evolves_boundary_variables or
147 : evolution::dg::system_has_boundary_variables_v<System>,
148 : "This boundary condition opts into evolving boundary variables, but "
149 : "the system's variables_tag has no ::Tags::BoundaryVariables entry.");
150 : // DemandOutgoingCharSpeeds should not need boundary evolved fields
151 : // as it does not write to external state.
152 : static_assert(
153 : not(bc_evolves_boundary_variables and
154 : BoundaryCondition::bc_type ==
155 : evolution::BoundaryConditions::Type::DemandOutgoingCharSpeeds),
156 : "A boundary condition that evolves boundary variables must not have "
157 : "`bc_type == DemandOutgoingCharSpeeds`.");
158 : using boundary_variables_tag_or_empty =
159 : tmpl::conditional_t<bc_evolves_boundary_variables,
160 : evolution::dg::boundary_variables_tag<System>,
161 : ::Tags::Variables<tmpl::list<>>>;
162 : using boundary_field_tags =
163 : typename boundary_variables_tag_or_empty::tags_list;
164 :
165 : // The interior fields `boundary_field_time_derivatives` need as input.
166 : // `boundary_field_time_derivatives` runs only on the physical pass.
167 : using bdry_evolved_field_projected_evolved_vars_tags = tmpl::conditional_t<
168 : bc_evolves_boundary_variables and not ComputeAuxiliary,
169 : get_boundary_field_time_derivatives_evolved_variables_tags_or_default_t<
170 : BoundaryCondition, tmpl::list<>>,
171 : tmpl::list<>>;
172 : using bdry_evolved_field_projected_primitive_tags = tmpl::conditional_t<
173 : bc_evolves_boundary_variables and not ComputeAuxiliary,
174 : get_boundary_field_time_derivatives_primitive_tags_or_default_t<
175 : BoundaryCondition, tmpl::list<>>,
176 : tmpl::list<>>;
177 : using bdry_evolved_field_projected_temporary_tags = tmpl::conditional_t<
178 : bc_evolves_boundary_variables and not ComputeAuxiliary,
179 : get_boundary_field_time_derivatives_temporary_tags_or_default_t<
180 : BoundaryCondition, tmpl::list<>>,
181 : tmpl::list<>>;
182 : static_assert(
183 : System::has_primitive_and_conservative_vars or
184 : tmpl::size<bdry_evolved_field_projected_primitive_tags>::value == 0,
185 : "A boundary condition cannot request "
186 : "`boundary_field_time_derivatives_primitive_tags` when the system has no "
187 : "primitive variables.");
188 :
189 : const Mesh<Dim - 1> face_mesh = volume_mesh.slice_away(direction.dimension());
190 : const size_t number_of_points_on_face = face_mesh.number_of_grid_points();
191 :
192 : // We figure out all the tags we need to project from the interior, both for
193 : // the boundary condition computation and for the boundary correction. We do
194 : // this by:
195 : // 1. get all interior tags for the boundary condition
196 : // 2. get all interior tags for the boundary correction (if ghost condition)
197 : // 3. combine these lists
198 : // 4. project from the interior
199 : //
200 : // Note: we only need to consider the boundary correction tags if a ghost
201 : // boundary condition is imposed.
202 :
203 : constexpr bool uses_ghost_condition =
204 : BoundaryCondition::bc_type ==
205 : evolution::BoundaryConditions::Type::Ghost or
206 : BoundaryCondition::bc_type ==
207 : evolution::BoundaryConditions::Type::GhostAndTimeDerivative;
208 : constexpr bool uses_time_derivative_condition =
209 : BoundaryCondition::bc_type ==
210 : evolution::BoundaryConditions::Type::TimeDerivative or
211 : BoundaryCondition::bc_type ==
212 : evolution::BoundaryConditions::Type::GhostAndTimeDerivative;
213 : // List that holds the inverse spatial metric if it's needed
214 : using inverse_spatial_metric_list =
215 : detail::inverse_spatial_metric_tag<System>;
216 : constexpr bool has_inv_spatial_metric =
217 : detail::has_inverse_spatial_metric_tag_v<System>;
218 :
219 : // Set up tags for boundary conditions
220 : using bcondition_interior_temp_tags =
221 : typename BoundaryCondition::dg_interior_temporary_tags;
222 : using bcondition_interior_prim_tags =
223 : detail::boundary_condition_primitive_tags<
224 : System::has_primitive_and_conservative_vars, BoundaryCondition>;
225 : using bcondition_interior_evolved_vars_tags =
226 : typename BoundaryCondition::dg_interior_evolved_variables_tags;
227 : using bcondition_interior_dt_evolved_vars_tags =
228 : detail::get_dt_vars_from_boundary_condition<BoundaryCondition>;
229 : using bcondition_interior_deriv_evolved_vars_tags =
230 : detail::get_deriv_vars_from_boundary_condition<BoundaryCondition>;
231 : using bcondition_interior_tags = tmpl::append<
232 : tmpl::conditional_t<has_inv_spatial_metric,
233 : tmpl::list<detail::NormalVector<Dim>>, tmpl::list<>>,
234 : bcondition_interior_evolved_vars_tags, bcondition_interior_prim_tags,
235 : bcondition_interior_temp_tags, bcondition_interior_dt_evolved_vars_tags,
236 : bcondition_interior_deriv_evolved_vars_tags>;
237 :
238 : // Set up tags for boundary correction
239 : using correction_temp_tags = tmpl::conditional_t<
240 : uses_ghost_condition,
241 : typename BoundaryCorrection::dg_package_data_temporary_tags,
242 : tmpl::list<>>;
243 : using correction_prim_tags = tmpl::conditional_t<
244 : uses_ghost_condition,
245 : detail::boundary_correction_primitive_tags<
246 : System::has_primitive_and_conservative_vars, BoundaryCorrection>,
247 : tmpl::list<>>;
248 : using correction_evolved_and_auxiliary_vars_tags =
249 : tmpl::conditional_t<uses_ghost_condition,
250 : tmpl::append<variables_tags, auxiliary_variables>,
251 : tmpl::list<>>;
252 :
253 : // Now combine the tags lists for each type of tag. These are all the tags
254 : // we need to project from the interior, excluding the inverse spatial
255 : // metric. They are the input to `dg_package_data` in the boundary
256 : // correction.
257 : using interior_temp_tags = tmpl::remove_duplicates<
258 : tmpl::append<bcondition_interior_temp_tags, correction_temp_tags,
259 : bdry_evolved_field_projected_temporary_tags>>;
260 : using interior_prim_tags = tmpl::remove_duplicates<
261 : tmpl::append<bcondition_interior_prim_tags, correction_prim_tags,
262 : bdry_evolved_field_projected_primitive_tags>>;
263 : using interior_vars_tags = tmpl::remove_duplicates<
264 : tmpl::append<correction_evolved_and_auxiliary_vars_tags,
265 : bcondition_interior_evolved_vars_tags,
266 : bdry_evolved_field_projected_evolved_vars_tags>>;
267 :
268 : constexpr bool needs_coordinates =
269 : tmpl::list_contains_v<interior_temp_tags,
270 : ::domain::Tags::Coordinates<Dim, Frame::Inertial>>;
271 :
272 : // List tags on the interior of the face. We list the exterior side
273 : // separately in the `else` branch of the if-constexpr where we actually use
274 : // the exterior fields.
275 : using fluxes_tags =
276 : tmpl::conditional_t<uses_ghost_condition,
277 : db::wrap_tags_in<::Tags::Flux, flux_variables,
278 : tmpl::size_t<Dim>, Frame::Inertial>,
279 : tmpl::list<>>;
280 : using tags_on_interior_face = tmpl::remove_duplicates<tmpl::append<
281 : fluxes_tags, interior_temp_tags, interior_prim_tags, interior_vars_tags,
282 : bcondition_interior_dt_evolved_vars_tags,
283 : bcondition_interior_deriv_evolved_vars_tags, inverse_spatial_metric_list,
284 : tmpl::list<detail::OneOverNormalVectorMagnitude,
285 : detail::NormalVector<Dim>>>>;
286 :
287 : Variables<tags_on_interior_face> interior_face_fields{
288 : number_of_points_on_face};
289 :
290 : // Perform projection into `interior_face_fields`. This also covers all the
291 : // fields for the exterior except for the time derivatives that might be
292 : // needed for Bjorhus/TimeDerivative boundary conditions.
293 : //
294 : // Note on the ordering of the data to project: if we are using a ghost
295 : // boundary condition with a boundary correction, then we know that all the
296 : // evolved variables are needed, whereas when using DemandOutgoingCharSpeeds
297 : // or Bjorhus boundary conditions none of the evolved variables might be
298 : // needed (or only some subset). Also, the way the typelist is assembled, the
299 : // evolved vars are guaranteed to be contiguous, but only if we are doing a
300 : // ghost boundary condition.
301 : if constexpr (uses_ghost_condition) {
302 : ::dg::project_contiguous_data_to_boundary(
303 : make_not_null(&interior_face_fields), volume_evolved_vars, volume_mesh,
304 : direction);
305 : if constexpr (tmpl::size<auxiliary_variables>::value != 0) {
306 : ASSERT(volume_auxiliary_variables != nullptr,
307 : "The auxiliary variables must be provided when the system has "
308 : "auxiliary variables.");
309 : ::dg::project_tensors_to_boundary<auxiliary_variables>(
310 : make_not_null(&interior_face_fields), *volume_auxiliary_variables,
311 : volume_mesh, direction);
312 : }
313 : } else {
314 : ::dg::project_tensors_to_boundary<interior_vars_tags>(
315 : make_not_null(&interior_face_fields), volume_evolved_vars, volume_mesh,
316 : direction);
317 : }
318 : if constexpr (tmpl::size<fluxes_tags>::value != 0) {
319 : ::dg::project_contiguous_data_to_boundary(
320 : make_not_null(&interior_face_fields), volume_fluxes, volume_mesh,
321 : direction);
322 : } else {
323 : (void)volume_fluxes;
324 : }
325 : using temp_tags_no_coordinates =
326 : tmpl::remove<interior_temp_tags,
327 : domain::Tags::Coordinates<Dim, Frame::Inertial>>;
328 : if constexpr (tmpl::size<tmpl::append<
329 : temp_tags_no_coordinates,
330 : detail::inverse_spatial_metric_tag<System>>>::value != 0) {
331 : ::dg::project_tensors_to_boundary<tmpl::append<
332 : temp_tags_no_coordinates, detail::inverse_spatial_metric_tag<System>>>(
333 : make_not_null(&interior_face_fields), volume_temporaries, volume_mesh,
334 : direction);
335 : }
336 : if constexpr (System::has_primitive_and_conservative_vars and
337 : tmpl::size<interior_prim_tags>::value != 0) {
338 : ASSERT(volume_primitive_variables != nullptr,
339 : "The volume primitive variables are not set even though the "
340 : "system has primitive variables.");
341 : ::dg::project_tensors_to_boundary<interior_prim_tags>(
342 : make_not_null(&interior_face_fields), *volume_primitive_variables,
343 : volume_mesh, direction);
344 : } else {
345 : (void)volume_primitive_variables;
346 : }
347 : if constexpr (tmpl::size<
348 : bcondition_interior_deriv_evolved_vars_tags>::value != 0) {
349 : ::dg::project_tensors_to_boundary<
350 : bcondition_interior_deriv_evolved_vars_tags>(
351 : make_not_null(&interior_face_fields), partial_derivs, volume_mesh,
352 : direction);
353 : }
354 : if constexpr (tmpl::size<bcondition_interior_dt_evolved_vars_tags>::value !=
355 : 0) {
356 : ::dg::project_tensors_to_boundary<bcondition_interior_dt_evolved_vars_tags>(
357 : make_not_null(&interior_face_fields), db::get<dt_variables_tag>(*box),
358 : volume_mesh, direction);
359 : }
360 :
361 : std::optional<tnsr::I<DataVector, Dim>> face_mesh_velocity{};
362 : if (volume_mesh_velocity.has_value()) {
363 : face_mesh_velocity = tnsr::I<DataVector, Dim>{number_of_points_on_face};
364 : ::dg::project_tensor_to_boundary(make_not_null(&*face_mesh_velocity),
365 : *volume_mesh_velocity, volume_mesh,
366 : direction);
367 : }
368 :
369 : // Normalize the normal vectors. We cache the unit normal covector For
370 : // flat geometry and static meshes.
371 : const auto normalize_normal_vectors =
372 : [&direction, mesh_is_moving = not moving_mesh_map.is_identity(),
373 : number_of_points_on_face, &volume_inverse_jacobian,
374 : &volume_mesh](const auto normal_covector_magnitude_in_direction_ptr,
375 : auto fields_on_face_ptr) {
376 : if (auto& normal_covector_quantity =
377 : *normal_covector_magnitude_in_direction_ptr;
378 : has_inv_spatial_metric or mesh_is_moving or
379 : not normal_covector_quantity.has_value()) {
380 : if (not normal_covector_quantity.has_value()) {
381 : normal_covector_quantity =
382 : Variables<tmpl::list<evolution::dg::Tags::MagnitudeOfNormal,
383 : evolution::dg::Tags::NormalCovector<Dim>>>{
384 : number_of_points_on_face};
385 : }
386 : tnsr::i<DataVector, Dim> volume_unnormalized_normal_covector{};
387 :
388 : for (size_t inertial_index = 0; inertial_index < Dim;
389 : ++inertial_index) {
390 : volume_unnormalized_normal_covector.get(inertial_index)
391 : .set_data_ref(
392 : const_cast<double*>( // NOLINT
393 : volume_inverse_jacobian
394 : .get(direction.dimension(), inertial_index)
395 : .data()),
396 : volume_mesh.number_of_grid_points());
397 : }
398 : ::dg::project_tensor_to_boundary(
399 : make_not_null(&get<evolution::dg::Tags::NormalCovector<Dim>>(
400 : *normal_covector_quantity)),
401 : volume_unnormalized_normal_covector, volume_mesh, direction);
402 :
403 : if (const double sign = direction.sign(); sign != 1.0) {
404 : for (auto& normal_covector_component :
405 : get<evolution::dg::Tags::NormalCovector<Dim>>(
406 : *normal_covector_quantity)) {
407 : normal_covector_component *= sign;
408 : }
409 : }
410 :
411 : detail::unit_normal_vector_and_covector_and_magnitude_impl<System>(
412 : make_not_null(&get<evolution::dg::Tags::MagnitudeOfNormal>(
413 : *normal_covector_quantity)),
414 : make_not_null(&get<evolution::dg::Tags::NormalCovector<Dim>>(
415 : *normal_covector_quantity)),
416 : fields_on_face_ptr,
417 : get<evolution::dg::Tags::NormalCovector<Dim>>(
418 : *normal_covector_quantity));
419 : }
420 : };
421 : // Normalize the outward facing normal vector on the interior side
422 : db::mutate<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>>(
423 : [&direction, &interior_face_fields, &normalize_normal_vectors](
424 : const auto normal_covector_and_magnitude_ptr) {
425 : normalize_normal_vectors(
426 : make_not_null(&normal_covector_and_magnitude_ptr->at(direction)),
427 : make_not_null(&interior_face_fields));
428 : },
429 : box);
430 :
431 : const tnsr::i<DataVector, Dim, Frame::Inertial>& interior_normal_covector =
432 : get<evolution::dg::Tags::NormalCovector<Dim>>(
433 : *db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>>(*box)
434 : .at(direction));
435 :
436 : if constexpr (needs_coordinates) {
437 : // Compute the coordinates on the interface
438 : get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(interior_face_fields) =
439 : moving_mesh_map(logical_to_grid_map(interface_logical_coordinates(
440 : face_mesh, direction)),
441 : time, functions_of_time);
442 : }
443 :
444 : if constexpr (BoundaryCondition::bc_type ==
445 : evolution::BoundaryConditions::Type::DemandOutgoingCharSpeeds) {
446 : // DemandOutgoingCharSpeeds boundary conditions only check that all
447 : // characteristic speeds are directed out of the element. If there are any
448 : // inward directed fields then the boundary condition should error.
449 : const auto apply_bc =
450 : [&boundary_condition, &face_mesh_velocity,
451 : &interior_normal_covector](const auto&... face_and_volume_args) {
452 : return boundary_condition.dg_demand_outgoing_char_speeds(
453 : face_mesh_velocity, interior_normal_covector,
454 : face_and_volume_args...);
455 : };
456 : const std::optional<std::string> error_message =
457 : apply_boundary_condition_impl(
458 : apply_bc, interior_face_fields, bcondition_interior_tags{},
459 : db::get<BoundaryConditionVolumeTags>(*box)...);
460 : if (error_message.has_value()) {
461 : ERROR(*error_message << "\n\nIn element:" << element.id()
462 : << "\nIn direction: " << direction);
463 : }
464 : return;
465 : }
466 :
467 : // We add the time derivative boundary conditions and lift the ghost boundary
468 : // conditions after both have been computed in case either depends on the
469 : // time derivatives in the volume projected on to the face.
470 :
471 : Variables<dt_variables_tags> dt_time_derivative_correction{};
472 : if constexpr (uses_time_derivative_condition and not ComputeAuxiliary) {
473 : dt_time_derivative_correction.initialize(number_of_points_on_face);
474 : auto apply_bc = [&boundary_condition, &dt_time_derivative_correction,
475 : &face_mesh_velocity, &interior_normal_covector](
476 : const auto&... interior_face_and_volume_args) {
477 : return boundary_condition.dg_time_derivative(
478 : make_not_null(&get<::Tags::dt<EvolvedVariablesTags>>(
479 : dt_time_derivative_correction))...,
480 : face_mesh_velocity, interior_normal_covector,
481 : interior_face_and_volume_args...);
482 : };
483 : const std::optional<std::string> error_message =
484 : apply_boundary_condition_impl(
485 : apply_bc, interior_face_fields, bcondition_interior_tags{},
486 : db::get<BoundaryConditionVolumeTags>(*box)...);
487 : if (error_message.has_value()) {
488 : ERROR(*error_message << "\n\nIn element:" << element.id()
489 : << "\nIn direction: " << direction);
490 : }
491 : } else {
492 : (void)dt_time_derivative_correction;
493 : }
494 :
495 : [[maybe_unused]] const Variables<boundary_field_tags>*
496 : stored_boundary_variables = nullptr;
497 : if constexpr (bc_evolves_boundary_variables) {
498 : stored_boundary_variables =
499 : &db::get<boundary_variables_tag_or_empty>(*box).variables().at(
500 : direction);
501 : }
502 : if constexpr (bc_evolves_boundary_variables and not ComputeAuxiliary) {
503 : using bdry_evolved_field_interior_tags =
504 : evolution::dg::boundary_field_time_derivatives_interior_tags<
505 : BoundaryCondition>;
506 : using dt_boundary_variables_tag =
507 : db::add_tag_prefix<::Tags::dt, boundary_variables_tag_or_empty>;
508 :
509 : db::mutate_apply<tmpl::list<dt_boundary_variables_tag>,
510 : tmpl::list<BoundaryConditionVolumeTags...>>(
511 : [&](const auto dt_boundary_variables, const auto&... volume_args) {
512 : auto& entry = dt_boundary_variables->variables().at(direction);
513 : ASSERT(
514 : entry.number_of_grid_points() == number_of_points_on_face,
515 : "The dt storage of the boundary-evolved variables for direction "
516 : << direction << " has " << entry.number_of_grid_points()
517 : << " grid points, but the face mesh has "
518 : << number_of_points_on_face
519 : << ". The entry must be sized to the face mesh.");
520 : const std::optional<std::string> error_message = tmpl::as_pack<
521 : boundary_field_tags>(
522 : [&]<typename... BoundaryFieldTags>(
523 : tmpl::type_<BoundaryFieldTags>... /*meta*/) {
524 : auto apply_boundary_field_time_derivatives =
525 : [&](const auto&... interior_face_and_volume_args) {
526 : return boundary_condition.boundary_field_time_derivatives(
527 : make_not_null(
528 : &get<::Tags::dt<BoundaryFieldTags>>(entry))...,
529 : face_mesh_velocity, interior_normal_covector,
530 : get<BoundaryFieldTags>(*stored_boundary_variables)...,
531 : interior_face_and_volume_args...);
532 : };
533 : return apply_boundary_condition_impl(
534 : apply_boundary_field_time_derivatives, interior_face_fields,
535 : bdry_evolved_field_interior_tags{}, volume_args...);
536 : });
537 : if (error_message.has_value()) {
538 : ERROR(*error_message << "\n\nIn element:" << element.id()
539 : << "\nIn direction: " << direction);
540 : }
541 : },
542 : box);
543 : }
544 :
545 : // Now we populate the fields on the exterior side of the face using the
546 : // boundary condition.
547 : // `auxiliary_variables` is included unconditionally: `dg_ghost` supplies the
548 : // exterior values of every field it declares in both the physical and the
549 : // auxiliary pass.
550 : using tags_on_exterior_face = tmpl::remove_duplicates<tmpl::append<
551 : variables_tags, auxiliary_variables, fluxes_tags, correction_temp_tags,
552 : correction_prim_tags, inverse_spatial_metric_list,
553 : tmpl::list<detail::OneOverNormalVectorMagnitude,
554 : detail::NormalVector<Dim>,
555 : evolution::dg::Tags::NormalCovector<Dim>>>>;
556 : Variables<tags_on_exterior_face> exterior_face_fields{
557 : number_of_points_on_face};
558 :
559 : const bool has_collocation_points_on_side =
560 : volume_mesh.quadrature(direction.dimension()) ==
561 : Spectral::Quadrature::GaussLobatto or
562 : (volume_mesh.quadrature(direction.dimension()) ==
563 : Spectral::Quadrature::GaussRadauUpper and
564 : direction.side() == Side::Upper);
565 :
566 : if constexpr (uses_ghost_condition) {
567 : using mortar_tags_list = tmpl::list<PackageFieldTags...>;
568 : using dg_package_data_projected_tags =
569 : tmpl::append<variables_tags, projected_auxiliary_vars_tags, fluxes_tags,
570 : correction_temp_tags, correction_prim_tags>;
571 :
572 : Variables<mortar_tags_list> internal_packaged_data{
573 : number_of_points_on_face};
574 : if constexpr (ComputeAuxiliary) {
575 : const double max_abs_char_speed_on_face =
576 : detail::dg_auxiliary_package_data<System>(
577 : make_not_null(&internal_packaged_data), boundary_correction,
578 : interior_face_fields, interior_normal_covector,
579 : face_mesh_velocity, dg_package_data_projected_tags{},
580 : db::get<PackageDataVolumeTags>(*box)...);
581 : (void)max_abs_char_speed_on_face;
582 : } else {
583 : const double max_abs_char_speed_on_face = detail::dg_package_data<System>(
584 : make_not_null(&internal_packaged_data), boundary_correction,
585 : interior_face_fields, interior_normal_covector, face_mesh_velocity,
586 : dg_package_data_projected_tags{},
587 : db::get<PackageDataVolumeTags>(*box)...);
588 : (void)max_abs_char_speed_on_face;
589 : }
590 :
591 : // Notes:
592 : // - We pass the outward directed normal vector normalized using the
593 : // interior variables to the boundary condition. This is because the
594 : // boundary condition should only need the normal vector for computing
595 : // things like reflecting BCs where the normal component of an interior
596 : // quantity is reversed.
597 : // - If needed, the boundary condition returns the inverse spatial metric on
598 : // the exterior side, which is then used to normalize the normal vector on
599 : // the exterior side. We need the exterior normal vector for computing
600 : // flux terms. The inverse spatial metric on the exterior side can be
601 : // equal to the inverse spatial metric on the interior side. This would be
602 : // true when, e.g. imposing reflecting boundary conditions.
603 : // - In addition to the evolved variables and fluxes, the boundary condition
604 : // must compute the `dg_packaged_data_temporary_tags` and the primitive
605 : // tags that the boundary correction needs.
606 : // - For systems with constraint damping parameters, the constraint damping
607 : // parameters are just copied from the projected values from the interior.
608 : // - If a boundary condition has boundary evolved fields, those fields are
609 : // passed as input arguments to dg_ghost. Their time derivatives are
610 : // computed by the boundary condition's method
611 : // boundary_field_time_derivatives.
612 : const auto boundary_field_values_args = [&stored_boundary_variables]() {
613 : if constexpr (bc_evolves_boundary_variables) {
614 : return tmpl::as_pack<boundary_field_tags>(
615 : [&stored_boundary_variables]<typename... BoundaryValueTags>(
616 : tmpl::type_<BoundaryValueTags>... /*meta*/) {
617 : return std::forward_as_tuple(
618 : get<BoundaryValueTags>(*stored_boundary_variables)...);
619 : });
620 : } else {
621 : (void)stored_boundary_variables;
622 : return std::tuple<>{};
623 : }
624 : }();
625 : auto apply_bc = [&boundary_condition, &exterior_face_fields,
626 : &face_mesh_velocity, &interior_normal_covector,
627 : &boundary_field_values_args](
628 : const auto&... interior_face_and_volume_args) {
629 : if constexpr (has_inv_spatial_metric) {
630 : return std::apply(
631 : [&](const auto&... boundary_field_values) {
632 : return boundary_condition.dg_ghost(
633 : make_not_null(&get<BoundaryCorrectionPackagedDataInputTags>(
634 : exterior_face_fields))...,
635 : make_not_null(
636 : &get<tmpl::front<
637 : detail::inverse_spatial_metric_tag<System>>>(
638 : exterior_face_fields)),
639 : face_mesh_velocity, interior_normal_covector,
640 : boundary_field_values..., interior_face_and_volume_args...);
641 : },
642 : boundary_field_values_args);
643 : } else {
644 : return std::apply(
645 : [&](const auto&... boundary_field_values) {
646 : return boundary_condition.dg_ghost(
647 : make_not_null(&get<BoundaryCorrectionPackagedDataInputTags>(
648 : exterior_face_fields))...,
649 : face_mesh_velocity, interior_normal_covector,
650 : boundary_field_values..., interior_face_and_volume_args...);
651 : },
652 : boundary_field_values_args);
653 : }
654 : };
655 : const std::optional<std::string> error_message =
656 : apply_boundary_condition_impl(
657 : apply_bc, interior_face_fields, bcondition_interior_tags{},
658 : db::get<BoundaryConditionVolumeTags>(*box)...);
659 : if (error_message.has_value()) {
660 : ERROR(*error_message << "\n\nIn element:" << element.id()
661 : << "\nIn direction: " << direction);
662 : }
663 : // Subtract mesh velocity from the _exterior_ fluxes
664 : if (face_mesh_velocity.has_value()) {
665 : tmpl::for_each<flux_variables>(
666 : [&face_mesh_velocity, &exterior_face_fields](auto tag_v) {
667 : // Modify fluxes for moving mesh
668 : using var_tag = typename decltype(tag_v)::type;
669 : using flux_var_tag =
670 : db::add_tag_prefix<::Tags::Flux, var_tag, tmpl::size_t<Dim>,
671 : Frame::Inertial>;
672 : auto& flux_var = get<flux_var_tag>(exterior_face_fields);
673 : const auto& var = get<var_tag>(exterior_face_fields);
674 : const auto& mesh_velocity = *face_mesh_velocity;
675 : // Loop over all independent components of flux_var
676 : for (size_t flux_var_storage_index = 0;
677 : flux_var_storage_index < flux_var.size();
678 : ++flux_var_storage_index) {
679 : // Get the flux variable's tensor index, e.g. (i,j) for a F^i of
680 : // the spatial velocity (or some other spatial tensor).
681 : const auto flux_var_tensor_index =
682 : flux_var.get_tensor_index(flux_var_storage_index);
683 : // Remove the first index from the flux tensor index, gets back
684 : // (j)
685 : const auto var_tensor_index =
686 : all_but_specified_element_of(flux_var_tensor_index, 0);
687 : // Set flux_index to (i)
688 : const size_t flux_index = gsl::at(flux_var_tensor_index, 0);
689 :
690 : // We now need to index flux(i,j) -= u(j) * v_g(i)
691 : flux_var[flux_var_storage_index] -=
692 : var.get(var_tensor_index) * mesh_velocity.get(flux_index);
693 : }
694 : });
695 : }
696 : // Now that we have computed the inverse spatial metric on the exterior, we
697 : // can compute the normalized normal (co)vector on the exterior side. If
698 : // there is no inverse spatial metric, then we just copy from the interior
699 : // and reverse the sign.
700 : for (size_t i = 0; i < Dim; ++i) {
701 : get<evolution::dg::Tags::NormalCovector<Dim>>(exterior_face_fields)
702 : .get(i) = -interior_normal_covector.get(i);
703 : }
704 : if constexpr (has_inv_spatial_metric) {
705 : const tnsr::II<DataVector, Dim, Frame::Inertial>& inv_spatial_metric =
706 : get<tmpl::front<inverse_spatial_metric_list>>(exterior_face_fields);
707 : tnsr::i<DataVector, Dim, Frame::Inertial>& exterior_normal_covector =
708 : get<evolution::dg::Tags::NormalCovector<Dim>>(exterior_face_fields);
709 : tnsr::I<DataVector, Dim, Frame::Inertial>& exterior_normal_vector =
710 : get<detail::NormalVector<Dim>>(exterior_face_fields);
711 :
712 : // Since the spatial metric is different on the exterior side of the
713 : // interface, we need to normalize the direction-reversed interior normal
714 : // vector using the exterior inverse spatial metric.
715 : for (size_t i = 0; i < Dim; ++i) {
716 : exterior_normal_vector.get(i) =
717 : get<0>(exterior_normal_covector) * inv_spatial_metric.get(i, 0);
718 : for (size_t j = 1; j < Dim; ++j) {
719 : exterior_normal_vector.get(i) +=
720 : exterior_normal_covector.get(j) * inv_spatial_metric.get(i, j);
721 : }
722 : }
723 : // Use detail::OneOverNormalVectorMagnitude as a buffer for the
724 : // magnitude. We don't need one over the normal magnitude on the
725 : // exterior side since we aren't lifting there.
726 : Scalar<DataVector>& magnitude =
727 : get<detail::OneOverNormalVectorMagnitude>(exterior_face_fields);
728 : dot_product(make_not_null(&magnitude), exterior_normal_covector,
729 : exterior_normal_vector);
730 : get(magnitude) = sqrt(get(magnitude));
731 : for (size_t i = 0; i < Dim; ++i) {
732 : exterior_normal_covector.get(i) /= get(magnitude);
733 : exterior_normal_vector.get(i) /= get(magnitude);
734 : }
735 : }
736 :
737 : // Package the external-side data for the boundary correction
738 : Variables<mortar_tags_list> external_packaged_data{
739 : number_of_points_on_face};
740 : if constexpr (ComputeAuxiliary) {
741 : detail::dg_auxiliary_package_data<System>(
742 : make_not_null(&external_packaged_data), boundary_correction,
743 : exterior_face_fields,
744 : get<evolution::dg::Tags::NormalCovector<Dim>>(exterior_face_fields),
745 : face_mesh_velocity, dg_package_data_projected_tags{},
746 : db::get<PackageDataVolumeTags>(*box)...);
747 : } else {
748 : detail::dg_package_data<System>(
749 : make_not_null(&external_packaged_data), boundary_correction,
750 : exterior_face_fields,
751 : get<evolution::dg::Tags::NormalCovector<Dim>>(exterior_face_fields),
752 : face_mesh_velocity, dg_package_data_projected_tags{},
753 : db::get<PackageDataVolumeTags>(*box)...);
754 : }
755 :
756 : typename tag_to_update::type boundary_corrections_on_face{
757 : number_of_points_on_face};
758 :
759 : // Compute boundary correction
760 : if constexpr (ComputeAuxiliary) {
761 : // The auxiliary boundary terms write one correction per auxiliary
762 : // variable (the buffer's own tags), not per evolved variable.
763 : boundary_correction.dg_auxiliary_boundary_terms(
764 : make_not_null(
765 : &get<AuxiliaryCorrectionTags>(boundary_corrections_on_face))...,
766 : get<PackageFieldTags>(internal_packaged_data)...,
767 : get<PackageFieldTags>(external_packaged_data)..., dg_formulation,
768 : get<BoundaryTermsVolumeTags>(*box)...);
769 : } else {
770 : boundary_correction.dg_boundary_terms(
771 : make_not_null(&get<::Tags::dt<EvolvedVariablesTags>>(
772 : boundary_corrections_on_face))...,
773 : get<PackageFieldTags>(internal_packaged_data)...,
774 : get<PackageFieldTags>(external_packaged_data)..., dg_formulation,
775 : get<BoundaryTermsVolumeTags>(*box)...);
776 : }
777 :
778 : // Lift the boundary correction
779 : const auto& magnitude_of_interior_face_normal =
780 : get<evolution::dg::Tags::MagnitudeOfNormal>(
781 : *db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>>(*box)
782 : .at(direction));
783 : if (has_collocation_points_on_side) {
784 : // The lift_flux function lifts only on the slice, it does not add
785 : // the contribution to the volume.
786 : ::dg::lift_flux(make_not_null(&boundary_corrections_on_face),
787 : volume_mesh.extents(direction.dimension()),
788 : magnitude_of_interior_face_normal,
789 : volume_mesh.basis(direction.dimension()));
790 :
791 : // Add the flux contribution to the volume data
792 : db::mutate<tag_to_update>(
793 : [&direction, &boundary_corrections_on_face,
794 : &volume_mesh](const auto vars_ptr) {
795 : add_slice_to_data(
796 : vars_ptr, boundary_corrections_on_face, volume_mesh.extents(),
797 : direction.dimension(),
798 : index_to_slice_at(volume_mesh.extents(), direction));
799 : },
800 : box);
801 : } else {
802 : // We are using Gauss points.
803 : //
804 : // Optimization note: eliminate allocations for volume and face det
805 : // jacobian. Should probably compute face det inv jacobian, then divide
806 : // (fewer grid points => fewer FLOPs).
807 : const DataVector volume_det_jacobian = 1.0 / get(volume_det_inv_jacobian);
808 :
809 : // Project the determinant of the Jacobian to the face. This could
810 : // be optimized by caching in the time-independent case.
811 : Scalar<DataVector> face_det_jacobian{face_mesh.number_of_grid_points()};
812 : const Matrix identity{};
813 : auto interpolation_matrices = make_array<Dim>(std::cref(identity));
814 : const std::pair<Matrix, Matrix>& matrices =
815 : Spectral::boundary_interpolation_matrices(
816 : volume_mesh.slice_through(direction.dimension()));
817 : gsl::at(interpolation_matrices, direction.dimension()) =
818 : direction.side() == Side::Upper ? matrices.second : matrices.first;
819 : apply_matrices(make_not_null(&get(face_det_jacobian)),
820 : interpolation_matrices, volume_det_jacobian,
821 : volume_mesh.extents());
822 :
823 : db::mutate<tag_to_update>(
824 : [&direction, &boundary_corrections_on_face, &face_det_jacobian,
825 : &magnitude_of_interior_face_normal, &volume_det_inv_jacobian,
826 : &volume_mesh](const auto vars_ptr) {
827 : ::dg::lift_boundary_terms_gauss_points(
828 : vars_ptr, volume_det_inv_jacobian, volume_mesh, direction,
829 : boundary_corrections_on_face, magnitude_of_interior_face_normal,
830 : face_det_jacobian);
831 : },
832 : box);
833 : }
834 : }
835 : // Add TimeDerivative correction to volume time derivatives.
836 : if constexpr (uses_time_derivative_condition and not ComputeAuxiliary) {
837 : if (has_collocation_points_on_side) {
838 : db::mutate<dt_variables_tag>(
839 : [&direction, &dt_time_derivative_correction,
840 : &volume_mesh](const auto dt_variables_ptr) {
841 : add_slice_to_data(
842 : dt_variables_ptr, dt_time_derivative_correction,
843 : volume_mesh.extents(), direction.dimension(),
844 : index_to_slice_at(volume_mesh.extents(), direction));
845 : },
846 : box);
847 : } else {
848 : db::mutate<dt_variables_tag>(
849 : [&direction, &dt_time_derivative_correction,
850 : &volume_mesh](const auto dt_variables_ptr) {
851 : ::dg::interpolate_dt_terms_gauss_points(
852 : dt_variables_ptr, volume_mesh, direction,
853 : dt_time_derivative_correction);
854 : },
855 : box);
856 : }
857 : }
858 : }
859 :
860 : /*!
861 : * \brief Applies the boundary conditions using the `boundary_correction`
862 : * on all external faces.
863 : *
864 : * A `tmpl::for_each` loop along with a `typeid` comparison checks which of the
865 : * known boundary conditions is being used. Since each direction can have a
866 : * different boundary condition, we must check each boundary condition in
867 : * each external direction.
868 : *
869 : * When `ComputeAuxiliary` is `true` the LDG auxiliary pass is applied instead
870 : * of the physical pass: ghost faces package data with the boundary
871 : * correction's `dg_auxiliary_*` interface and lift the resulting correction
872 : * into the auxiliary-variable storage `::Tags::Variables<auxiliary_variables>`
873 : * rather than into the time derivatives, and time-derivative boundary
874 : * conditions are skipped.
875 : */
876 : template <typename System, size_t Dim, typename VariablesTag,
877 : bool ComputeAuxiliary = false, typename DbTagsList,
878 : typename BoundaryCorrection>
879 : void apply_boundary_conditions_on_all_external_faces(
880 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
881 : const BoundaryCorrection& boundary_correction,
882 : const Variables<
883 : typename System::compute_volume_time_derivative_terms::temporary_tags>&
884 : temporaries,
885 : const Variables<
886 : db::wrap_tags_in<::Tags::Flux, typename System::flux_variables,
887 : tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes,
888 : const Variables<
889 : db::wrap_tags_in<::Tags::deriv, typename System::gradient_variables,
890 : tmpl::size_t<Dim>, Frame::Inertial>>& partial_derivs,
891 : const Variables<detail::get_primitive_vars_tags_from_system<System>>* const
892 : primitive_vars) {
893 : using factory_classes =
894 : typename std::decay_t<decltype(db::get<Parallel::Tags::Metavariables>(
895 : *box))>::factory_creation::factory_classes;
896 :
897 : using derived_boundary_conditions = tmpl::remove_if<
898 : tmpl::at<factory_classes, typename System::boundary_conditions_base>,
899 : tmpl::or_<
900 : std::is_base_of<domain::BoundaryConditions::MarkAsCartoon, tmpl::_1>,
901 : std::is_base_of<domain::BoundaryConditions::MarkAsNone, tmpl::_1>,
902 : std::is_base_of<domain::BoundaryConditions::MarkAsPeriodic,
903 : tmpl::_1>>>;
904 :
905 : using variables_tag = VariablesTag;
906 : using flux_variables = typename System::flux_variables;
907 : using fluxes_tags = db::wrap_tags_in<::Tags::Flux, flux_variables,
908 : tmpl::size_t<Dim>, Frame::Inertial>;
909 : using auxiliary_variables =
910 : get_auxiliary_variables_or_default_t<System, tmpl::list<>>;
911 :
912 : const Element<Dim>& element = db::get<domain::Tags::Element<Dim>>(*box);
913 : size_t number_of_boundaries_left = element.external_boundaries().size();
914 :
915 : if (number_of_boundaries_left == 0) {
916 : return;
917 : }
918 :
919 : const Variables<auxiliary_variables>* volume_auxiliary_variables = nullptr;
920 : if constexpr (tmpl::size<auxiliary_variables>::value != 0) {
921 : volume_auxiliary_variables =
922 : &db::get<::Tags::Variables<auxiliary_variables>>(*box);
923 : }
924 :
925 : const auto& external_boundary_conditions =
926 : db::get<domain::Tags::ExternalBoundaryConditions<Dim>>(*box).at(
927 : element.id().block_id());
928 :
929 : // Error if any cartoon-type BC is used with an incompatible mesh.
930 : // Cartoon BCs are coupled only to ZernikeB1 basis elements, which do not need
931 : // a DG boundary condition but do require one with subcell
932 : if constexpr (domain::BoundaryConditions::detail::has_cartoon_bc_v<
933 : tmpl::at<factory_classes,
934 : typename System::boundary_conditions_base>>) {
935 : const auto& mesh = db::get<::domain::Tags::Mesh<Dim>>(*box);
936 : if (not domain::BoundaryConditions::dg_mesh_is_cartoon_compatible(mesh)) {
937 : for (const Direction<Dim>& direction : element.external_boundaries()) {
938 : if (domain::BoundaryConditions::is_cartoon(
939 : external_boundary_conditions.at(direction))) {
940 : ERROR(
941 : "You might have used a Cartoon boundary condition on an external "
942 : "boundary condition. Alternatively and less likely, there is a "
943 : "bug. The problematic BC is in direction "
944 : << direction << ", the mesh is: " << mesh);
945 : }
946 : }
947 : }
948 : }
949 :
950 : tmpl::for_each<
951 : derived_boundary_conditions>([&boundary_correction, &box, &element,
952 : &external_boundary_conditions,
953 : &number_of_boundaries_left, &partial_derivs,
954 : &primitive_vars, &temporaries,
955 : &volume_fluxes, volume_auxiliary_variables](
956 : auto derived_boundary_condition_v) {
957 : using DerivedBoundaryCondition =
958 : tmpl::type_from<decltype(derived_boundary_condition_v)>;
959 :
960 : if (number_of_boundaries_left == 0) {
961 : return;
962 : }
963 :
964 : for (const Direction<Dim>& direction : element.external_boundaries()) {
965 : const auto& boundary_condition =
966 : *external_boundary_conditions.at(direction);
967 : if (typeid(boundary_condition) == typeid(DerivedBoundaryCondition)) {
968 : // Select auxiliary or physical tags based on ComputeAuxiliary
969 : using package_data_volume_tags = tmpl::conditional_t<
970 : ComputeAuxiliary,
971 : get_dg_auxiliary_package_data_volume_tags_or_default_t<
972 : BoundaryCorrection, tmpl::list<>>,
973 : typename BoundaryCorrection::dg_package_data_volume_tags>;
974 : using package_field_tags = tmpl::conditional_t<
975 : ComputeAuxiliary,
976 : get_dg_auxiliary_package_field_tags_or_default_t<BoundaryCorrection,
977 : tmpl::list<>>,
978 : typename BoundaryCorrection::dg_package_field_tags>;
979 : using boundary_terms_volume_tags = tmpl::conditional_t<
980 : ComputeAuxiliary,
981 : get_dg_auxiliary_boundary_terms_volume_tags_or_default_t<
982 : BoundaryCorrection, tmpl::list<>>,
983 : typename BoundaryCorrection::dg_boundary_terms_volume_tags>;
984 : using package_data_temp_tags = tmpl::conditional_t<
985 : ComputeAuxiliary,
986 : get_dg_auxiliary_package_data_temporary_tags_or_default_t<
987 : BoundaryCorrection, tmpl::list<>>,
988 : typename BoundaryCorrection::dg_package_data_temporary_tags>;
989 : detail::apply_boundary_condition_on_face<System, Dim, ComputeAuxiliary>(
990 : box, boundary_correction,
991 : dynamic_cast<const DerivedBoundaryCondition&>(boundary_condition),
992 : direction, db::get<variables_tag>(*box), volume_auxiliary_variables,
993 : volume_fluxes, partial_derivs, temporaries, primitive_vars,
994 : db::get<::dg::Tags::Formulation>(*box),
995 : db::get<::domain::Tags::Mesh<Dim>>(*box),
996 : db::get<::domain::Tags::Element<Dim>>(*box),
997 : db::get<::domain::Tags::ElementMap<Dim, Frame::Grid>>(*box),
998 : db::get<::domain::CoordinateMaps::Tags::CoordinateMap<
999 : Dim, Frame::Grid, Frame::Inertial>>(*box),
1000 : db::get<::Tags::Time>(*box),
1001 : db::get<::domain::Tags::FunctionsOfTime>(*box),
1002 : db::get<::domain::Tags::MeshVelocity<Dim>>(*box),
1003 : db::get<::domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
1004 : Frame::Inertial>>(*box),
1005 : db::get<::domain::Tags::DetInvJacobian<Frame::ElementLogical,
1006 : Frame::Inertial>>(*box),
1007 : package_data_volume_tags{}, package_field_tags{},
1008 : boundary_terms_volume_tags{},
1009 : tmpl::remove_duplicates<tmpl::append<
1010 : typename variables_tag::tags_list, auxiliary_variables,
1011 : fluxes_tags, package_data_temp_tags,
1012 : typename detail::get_primitive_vars<
1013 : System::has_primitive_and_conservative_vars>::
1014 : template f<BoundaryCorrection>>>{},
1015 : typename DerivedBoundaryCondition::dg_gridless_tags{},
1016 : auxiliary_variables{});
1017 : --number_of_boundaries_left;
1018 : }
1019 : if (number_of_boundaries_left == 0) {
1020 : return;
1021 : }
1022 : }
1023 : });
1024 : }
1025 : } // namespace evolution::dg::Actions::detail
|
