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 "Domain/Block.hpp"
23 : #include "Domain/BoundaryConditions/None.hpp"
24 : #include "Domain/BoundaryConditions/Periodic.hpp"
25 : #include "Domain/CoordinateMaps/CoordinateMap.hpp"
26 : #include "Domain/Creators/Tags/ExternalBoundaryConditions.hpp"
27 : #include "Domain/Domain.hpp"
28 : #include "Domain/ElementMap.hpp"
29 : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
30 : #include "Domain/FunctionsOfTime/Tags.hpp"
31 : #include "Domain/InterfaceLogicalCoordinates.hpp"
32 : #include "Domain/Structure/Direction.hpp"
33 : #include "Domain/Structure/Element.hpp"
34 : #include "Domain/Structure/ElementId.hpp"
35 : #include "Domain/Tags.hpp"
36 : #include "Domain/TagsTimeDependent.hpp"
37 : #include "Evolution/BoundaryConditions/Type.hpp"
38 : #include "Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivativeHelpers.hpp"
39 : #include "Evolution/DiscontinuousGalerkin/Actions/NormalCovectorAndMagnitude.hpp"
40 : #include "Evolution/DiscontinuousGalerkin/Actions/PackageDataImpl.hpp"
41 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
42 : #include "NumericalAlgorithms/DiscontinuousGalerkin/InterpolateFromBoundary.hpp"
43 : #include "NumericalAlgorithms/DiscontinuousGalerkin/LiftFlux.hpp"
44 : #include "NumericalAlgorithms/DiscontinuousGalerkin/LiftFromBoundary.hpp"
45 : #include "NumericalAlgorithms/DiscontinuousGalerkin/ProjectToBoundary.hpp"
46 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Tags/Formulation.hpp"
47 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
48 : #include "NumericalAlgorithms/Spectral/BoundaryInterpolationMatrices.hpp"
49 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
50 : #include "NumericalAlgorithms/Spectral/Quadrature.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 : const bool has_collocation_points_on_side =
428 : volume_mesh.quadrature(direction.dimension()) ==
429 : Spectral::Quadrature::GaussLobatto or
430 : (volume_mesh.quadrature(direction.dimension()) ==
431 : Spectral::Quadrature::GaussRadauUpper and
432 : direction.side() == Side::Upper);
433 :
434 : if constexpr (uses_ghost_condition) {
435 : using mortar_tags_list = tmpl::list<PackageFieldTags...>;
436 : using dg_package_data_projected_tags =
437 : tmpl::append<variables_tags, fluxes_tags, correction_temp_tags,
438 : correction_prim_tags>;
439 :
440 : Variables<mortar_tags_list> internal_packaged_data{
441 : number_of_points_on_face};
442 : const double max_abs_char_speed_on_face = detail::dg_package_data<System>(
443 : make_not_null(&internal_packaged_data), boundary_correction,
444 : interior_face_fields, interior_normal_covector, face_mesh_velocity,
445 : dg_package_data_projected_tags{},
446 : db::get<PackageDataVolumeTags>(*box)...);
447 : (void)max_abs_char_speed_on_face;
448 :
449 : // Notes:
450 : // - we pass the outward directed normal vector normalized using the
451 : // interior variables to the boundary condition. This is because the
452 : // boundary condition should only need the normal vector for computing
453 : // things like reflecting BCs where the normal component of an interior
454 : // quantity is reversed.
455 : // - if needed, the boundary condition returns the inverse spatial metric on
456 : // the exterior side, which is then used to normalize the normal vector on
457 : // the exterior side. We need the exterior normal vector for computing
458 : // flux terms. The inverse spatial metric on the exterior side can be
459 : // equal to the inverse spatial metric on the interior side. This would be
460 : // true when, e.g. imposing reflecting boundary conditions.
461 : // - in addition to the evolved variables and fluxes, the boundary condition
462 : // must compute the `dg_packaged_data_temporary_tags` and the primitive
463 : // tags that the boundary correction needs.
464 : // - For systems with constraint damping parameters, the constraint damping
465 : // parameters are just copied from the projected values from the interior.
466 : auto apply_bc = [&boundary_condition, &exterior_face_fields,
467 : &face_mesh_velocity, &interior_normal_covector](
468 : const auto&... interior_face_and_volume_args) {
469 : if constexpr (has_inv_spatial_metric) {
470 : return boundary_condition.dg_ghost(
471 : make_not_null(&get<BoundaryCorrectionPackagedDataInputTags>(
472 : exterior_face_fields))...,
473 : make_not_null(
474 : &get<tmpl::front<detail::inverse_spatial_metric_tag<System>>>(
475 : exterior_face_fields)),
476 : face_mesh_velocity, interior_normal_covector,
477 : interior_face_and_volume_args...);
478 : } else {
479 : return boundary_condition.dg_ghost(
480 : make_not_null(&get<BoundaryCorrectionPackagedDataInputTags>(
481 : exterior_face_fields))...,
482 : face_mesh_velocity, interior_normal_covector,
483 : interior_face_and_volume_args...);
484 : }
485 : };
486 : const std::optional<std::string> error_message =
487 : apply_boundary_condition_impl(
488 : apply_bc, interior_face_fields, bcondition_interior_tags{},
489 : db::get<BoundaryConditionVolumeTags>(*box)...);
490 : if (error_message.has_value()) {
491 : ERROR(*error_message << "\n\nIn element:" << element.id()
492 : << "\nIn direction: " << direction);
493 : }
494 : // Subtract mesh velocity from the _exterior_ fluxes
495 : if (face_mesh_velocity.has_value()) {
496 : tmpl::for_each<flux_variables>(
497 : [&face_mesh_velocity, &exterior_face_fields](auto tag_v) {
498 : // Modify fluxes for moving mesh
499 : using var_tag = typename decltype(tag_v)::type;
500 : using flux_var_tag =
501 : db::add_tag_prefix<::Tags::Flux, var_tag, tmpl::size_t<Dim>,
502 : Frame::Inertial>;
503 : auto& flux_var = get<flux_var_tag>(exterior_face_fields);
504 : const auto& var = get<var_tag>(exterior_face_fields);
505 : const auto& mesh_velocity = *face_mesh_velocity;
506 : // Loop over all independent components of flux_var
507 : for (size_t flux_var_storage_index = 0;
508 : flux_var_storage_index < flux_var.size();
509 : ++flux_var_storage_index) {
510 : // Get the flux variable's tensor index, e.g. (i,j) for a F^i of
511 : // the spatial velocity (or some other spatial tensor).
512 : const auto flux_var_tensor_index =
513 : flux_var.get_tensor_index(flux_var_storage_index);
514 : // Remove the first index from the flux tensor index, gets back
515 : // (j)
516 : const auto var_tensor_index =
517 : all_but_specified_element_of(flux_var_tensor_index, 0);
518 : // Set flux_index to (i)
519 : const size_t flux_index = gsl::at(flux_var_tensor_index, 0);
520 :
521 : // We now need to index flux(i,j) -= u(j) * v_g(i)
522 : flux_var[flux_var_storage_index] -=
523 : var.get(var_tensor_index) * mesh_velocity.get(flux_index);
524 : }
525 : });
526 : }
527 : // Now that we have computed the inverse spatial metric on the exterior, we
528 : // can compute the normalized normal (co)vector on the exterior side. If
529 : // there is no inverse spatial metric, then we just copy from the interior
530 : // and reverse the sign.
531 : for (size_t i = 0; i < Dim; ++i) {
532 : get<evolution::dg::Tags::NormalCovector<Dim>>(exterior_face_fields)
533 : .get(i) = -interior_normal_covector.get(i);
534 : }
535 : if constexpr (has_inv_spatial_metric) {
536 : const tnsr::II<DataVector, Dim, Frame::Inertial>& inv_spatial_metric =
537 : get<tmpl::front<inverse_spatial_metric_list>>(exterior_face_fields);
538 : tnsr::i<DataVector, Dim, Frame::Inertial>& exterior_normal_covector =
539 : get<evolution::dg::Tags::NormalCovector<Dim>>(exterior_face_fields);
540 : tnsr::I<DataVector, Dim, Frame::Inertial>& exterior_normal_vector =
541 : get<detail::NormalVector<Dim>>(exterior_face_fields);
542 :
543 : // Since the spatial metric is different on the exterior side of the
544 : // interface, we need to normalize the direction-reversed interior normal
545 : // vector using the exterior inverse spatial metric.
546 : for (size_t i = 0; i < Dim; ++i) {
547 : exterior_normal_vector.get(i) =
548 : get<0>(exterior_normal_covector) * inv_spatial_metric.get(i, 0);
549 : for (size_t j = 1; j < Dim; ++j) {
550 : exterior_normal_vector.get(i) +=
551 : exterior_normal_covector.get(j) * inv_spatial_metric.get(i, j);
552 : }
553 : }
554 : // Use detail::OneOverNormalVectorMagnitude as a buffer for the
555 : // magnitude. We don't need one over the normal magnitude on the
556 : // exterior side since we aren't lifting there.
557 : Scalar<DataVector>& magnitude =
558 : get<detail::OneOverNormalVectorMagnitude>(exterior_face_fields);
559 : dot_product(make_not_null(&magnitude), exterior_normal_covector,
560 : exterior_normal_vector);
561 : get(magnitude) = sqrt(get(magnitude));
562 : for (size_t i = 0; i < Dim; ++i) {
563 : exterior_normal_covector.get(i) /= get(magnitude);
564 : exterior_normal_vector.get(i) /= get(magnitude);
565 : }
566 : }
567 :
568 : // Package the external-side data for the boundary correction
569 : Variables<mortar_tags_list> external_packaged_data{
570 : number_of_points_on_face};
571 : detail::dg_package_data<System>(
572 : make_not_null(&external_packaged_data), boundary_correction,
573 : exterior_face_fields,
574 : get<evolution::dg::Tags::NormalCovector<Dim>>(exterior_face_fields),
575 : face_mesh_velocity, dg_package_data_projected_tags{},
576 : db::get<PackageDataVolumeTags>(*box)...);
577 :
578 : Variables<dt_variables_tags> boundary_corrections_on_face{
579 : number_of_points_on_face};
580 :
581 : // Compute boundary correction
582 : boundary_correction.dg_boundary_terms(
583 : make_not_null(&get<::Tags::dt<EvolvedVariablesTags>>(
584 : boundary_corrections_on_face))...,
585 : get<PackageFieldTags>(internal_packaged_data)...,
586 : get<PackageFieldTags>(external_packaged_data)..., dg_formulation,
587 : get<BoundaryTermsVolumeTags>(*box)...);
588 :
589 : // Lift the boundary correction
590 : const auto& magnitude_of_interior_face_normal =
591 : get<evolution::dg::Tags::MagnitudeOfNormal>(
592 : *db::get<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>>(*box)
593 : .at(direction));
594 : if (has_collocation_points_on_side) {
595 : // The lift_flux function lifts only on the slice, it does not add
596 : // the contribution to the volume.
597 : ::dg::lift_flux(make_not_null(&boundary_corrections_on_face),
598 : volume_mesh.extents(direction.dimension()),
599 : magnitude_of_interior_face_normal,
600 : volume_mesh.basis(direction.dimension()));
601 :
602 : // Add the flux contribution to the volume data
603 : db::mutate<dt_variables_tag>(
604 : [&direction, &boundary_corrections_on_face,
605 : &volume_mesh](const auto dt_variables_ptr) {
606 : add_slice_to_data(
607 : dt_variables_ptr, boundary_corrections_on_face,
608 : volume_mesh.extents(), direction.dimension(),
609 : index_to_slice_at(volume_mesh.extents(), direction));
610 : },
611 : box);
612 : } else {
613 : // We are using Gauss points.
614 : //
615 : // Optimization note: eliminate allocations for volume and face det
616 : // jacobian. Should probably compute face det inv jacobian, then divide
617 : // (fewer grid points => fewer FLOPs).
618 : const DataVector volume_det_jacobian = 1.0 / get(volume_det_inv_jacobian);
619 :
620 : // Project the determinant of the Jacobian to the face. This could
621 : // be optimized by caching in the time-independent case.
622 : Scalar<DataVector> face_det_jacobian{face_mesh.number_of_grid_points()};
623 : const Matrix identity{};
624 : auto interpolation_matrices = make_array<Dim>(std::cref(identity));
625 : const std::pair<Matrix, Matrix>& matrices =
626 : Spectral::boundary_interpolation_matrices(
627 : volume_mesh.slice_through(direction.dimension()));
628 : gsl::at(interpolation_matrices, direction.dimension()) =
629 : direction.side() == Side::Upper ? matrices.second : matrices.first;
630 : apply_matrices(make_not_null(&get(face_det_jacobian)),
631 : interpolation_matrices, volume_det_jacobian,
632 : volume_mesh.extents());
633 :
634 : db::mutate<dt_variables_tag>(
635 : [&direction, &boundary_corrections_on_face, &face_det_jacobian,
636 : &magnitude_of_interior_face_normal, &volume_det_inv_jacobian,
637 : &volume_mesh](const auto dt_variables_ptr) {
638 : ::dg::lift_boundary_terms_gauss_points(
639 : dt_variables_ptr, volume_det_inv_jacobian, volume_mesh,
640 : direction, boundary_corrections_on_face,
641 : magnitude_of_interior_face_normal, face_det_jacobian);
642 : },
643 : box);
644 : }
645 : }
646 : // Add TimeDerivative correction to volume time derivatives.
647 : if constexpr (uses_time_derivative_condition) {
648 : if (has_collocation_points_on_side) {
649 : db::mutate<dt_variables_tag>(
650 : [&direction, &dt_time_derivative_correction,
651 : &volume_mesh](const auto dt_variables_ptr) {
652 : add_slice_to_data(
653 : dt_variables_ptr, dt_time_derivative_correction,
654 : volume_mesh.extents(), direction.dimension(),
655 : index_to_slice_at(volume_mesh.extents(), direction));
656 : },
657 : box);
658 : } else {
659 : db::mutate<dt_variables_tag>(
660 : [&direction, &dt_time_derivative_correction,
661 : &volume_mesh](const auto dt_variables_ptr) {
662 : ::dg::interpolate_dt_terms_gauss_points(
663 : dt_variables_ptr, volume_mesh, direction,
664 : dt_time_derivative_correction);
665 : },
666 : box);
667 : }
668 : }
669 : }
670 :
671 : /*!
672 : * \brief Applies the boundary conditions using the `boundary_correction`
673 : * on all external faces.
674 : *
675 : * A `tmpl::for_each` loop along with a `typeid` comparison checks which of the
676 : * known boundary conditions is being used. Since each direction can have a
677 : * different boundary condition, we must check each boundary condition in
678 : * each external direction.
679 : */
680 : template <typename System, size_t Dim, typename DbTagsList,
681 : typename BoundaryCorrection>
682 : void apply_boundary_conditions_on_all_external_faces(
683 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
684 : const BoundaryCorrection& boundary_correction,
685 : const Variables<
686 : typename System::compute_volume_time_derivative_terms::temporary_tags>&
687 : temporaries,
688 : const Variables<
689 : db::wrap_tags_in<::Tags::Flux, typename System::flux_variables,
690 : tmpl::size_t<Dim>, Frame::Inertial>>& volume_fluxes,
691 : const Variables<
692 : db::wrap_tags_in<::Tags::deriv, typename System::gradient_variables,
693 : tmpl::size_t<Dim>, Frame::Inertial>>& partial_derivs,
694 : const Variables<detail::get_primitive_vars_tags_from_system<System>>* const
695 : primitive_vars) {
696 : using factory_classes =
697 : typename std::decay_t<decltype(db::get<Parallel::Tags::Metavariables>(
698 : *box))>::factory_creation::factory_classes;
699 : using derived_boundary_conditions = tmpl::remove_if<
700 : tmpl::at<factory_classes, typename System::boundary_conditions_base>,
701 : tmpl::or_<
702 : std::is_base_of<domain::BoundaryConditions::MarkAsPeriodic, tmpl::_1>,
703 : std::is_base_of<domain::BoundaryConditions::MarkAsNone, tmpl::_1>>>;
704 :
705 : using variables_tag = typename System::variables_tag;
706 : using flux_variables = typename System::flux_variables;
707 : using fluxes_tags = db::wrap_tags_in<::Tags::Flux, flux_variables,
708 : tmpl::size_t<Dim>, Frame::Inertial>;
709 :
710 : const Element<Dim>& element = db::get<domain::Tags::Element<Dim>>(*box);
711 : size_t number_of_boundaries_left = element.external_boundaries().size();
712 :
713 : if (number_of_boundaries_left == 0) {
714 : return;
715 : }
716 :
717 : const auto& external_boundary_conditions =
718 : db::get<domain::Tags::ExternalBoundaryConditions<Dim>>(*box).at(
719 : element.id().block_id());
720 : tmpl::for_each<derived_boundary_conditions>(
721 : [&boundary_correction, &box, &element, &external_boundary_conditions,
722 : &number_of_boundaries_left, &partial_derivs, &primitive_vars,
723 : &temporaries, &volume_fluxes](auto derived_boundary_condition_v) {
724 : using DerivedBoundaryCondition =
725 : tmpl::type_from<decltype(derived_boundary_condition_v)>;
726 :
727 : if (number_of_boundaries_left == 0) {
728 : return;
729 : }
730 :
731 : for (const Direction<Dim>& direction : element.external_boundaries()) {
732 : const auto& boundary_condition =
733 : *external_boundary_conditions.at(direction);
734 : if (typeid(boundary_condition) == typeid(DerivedBoundaryCondition)) {
735 : detail::apply_boundary_condition_on_face<System>(
736 : box, boundary_correction,
737 : dynamic_cast<const DerivedBoundaryCondition&>(
738 : boundary_condition),
739 : direction, db::get<variables_tag>(*box), volume_fluxes,
740 : partial_derivs, temporaries, primitive_vars,
741 : db::get<::dg::Tags::Formulation>(*box),
742 : db::get<::domain::Tags::Mesh<Dim>>(*box),
743 : db::get<::domain::Tags::Element<Dim>>(*box),
744 : db::get<::domain::Tags::ElementMap<Dim, Frame::Grid>>(*box),
745 : db::get<::domain::CoordinateMaps::Tags::CoordinateMap<
746 : Dim, Frame::Grid, Frame::Inertial>>(*box),
747 : db::get<::Tags::Time>(*box),
748 : db::get<::domain::Tags::FunctionsOfTime>(*box),
749 : db::get<::domain::Tags::MeshVelocity<Dim>>(*box),
750 : db::get<::domain::Tags::InverseJacobian<
751 : Dim, Frame::ElementLogical, Frame::Inertial>>(*box),
752 : db::get<::domain::Tags::DetInvJacobian<Frame::ElementLogical,
753 : Frame::Inertial>>(*box),
754 : typename BoundaryCorrection::dg_package_data_volume_tags{},
755 : typename BoundaryCorrection::dg_package_field_tags{},
756 : typename BoundaryCorrection::dg_boundary_terms_volume_tags{},
757 : tmpl::remove_duplicates<tmpl::append<
758 : typename variables_tag::tags_list, fluxes_tags,
759 : typename BoundaryCorrection::dg_package_data_temporary_tags,
760 : typename detail::get_primitive_vars<
761 : System::has_primitive_and_conservative_vars>::
762 : template f<BoundaryCorrection>>>{},
763 : typename DerivedBoundaryCondition::dg_gridless_tags{});
764 : --number_of_boundaries_left;
765 : }
766 : if (number_of_boundaries_left == 0) {
767 : return;
768 : }
769 : }
770 : });
771 : }
772 : } // namespace evolution::dg::Actions::detail
|
