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