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