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