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