CorrectPackagedData.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <algorithm>
7 #include <boost/functional/hash.hpp>
8 #include <cstddef>
9 #include <unordered_map>
10 #include <utility>
11 
12 #include "DataStructures/DataVector.hpp"
13 #include "DataStructures/Index.hpp"
14 #include "DataStructures/SliceIterator.hpp"
20 #include "Evolution/DgSubcell/Projection.hpp"
21 #include "Evolution/DiscontinuousGalerkin/MortarData.hpp"
24 #include "Utilities/Gsl.hpp"
25 #include "Utilities/TMPL.hpp"
26 
27 namespace evolution::dg::subcell {
28 /*!
29  * \brief Project the DG package data to the subcells. Data received from a
30  * neighboring element doing DG is always projected, while the data we sent to
31  * our neighbors before doing a rollback from DG to subcell is only projected if
32  * `OverwriteInternalMortarData` is `true`.
33  *
34  * In order for the hybrid DG-FD/FV scheme to be conservative between elements
35  * using DG and elements using subcell, the boundary terms must be the same on
36  * both elements. In practice this means the boundary corrections \f$G+D\f$ must
37  * be computed on the same grid. Consider the element doing subcell which
38  * receives data from an element doing DG. In this case the DG element's
39  * ingredients going into \f$G+D\f$ are projected to the subcells and then
40  * \f$G+D\f$ are computed on the subcells. Similarly, for strict conservation
41  * the element doing DG must first project the data it sent to the neighbor to
42  * the subcells, then compute \f$G+D\f$ on the subcells, and finally reconstrct
43  * \f$G+D\f$ back to the DG grid before lifting \f$G+D\f$ to the volume.
44  *
45  * This function updates the `packaged_data` (ingredients into \f$G+D\f$)
46  * received by an element doing subcell by projecting the neighbor's DG data
47  * onto the subcells. Note that this is only half of what is required for strict
48  * conservation, the DG element must also compute \f$G+D\f$ on the subcells.
49  * Note that we currently do not perform the other half of the correction
50  * needed to be strictly conservative.
51  *
52  * If we are retaking a time step after the DG step failed then maintaining
53  * conservation requires additional care. If `OverwriteInternalMortarData` is
54  * `true` then the local (the element switching from DG to subcell) ingredients
55  * into \f$G+D\f$ are projected and overwrite the data computed from
56  * the FD reconstruction to the interface. However, even this is insufficient to
57  * guarantee conservation. To guarantee conservation (which we do not currently
58  * do) the correction \f$G+D\f$ must be computed on the DG grid and then
59  * projected to the subcells.
60  *
61  * Note that our practical experience shows that since the DG-subcell hybrid
62  * scheme switches to the subcell solver _before_ the local solution contains
63  * discontinuities, strict conservation is not necessary between DG and FD/FV
64  * regions. This was also observed with a block-adaptive finite difference AMR
65  * code \cite CHEN2016604
66  */
67 template <bool OverwriteInternalMortarData, size_t Dim,
68  typename DgPackageFieldTags>
70  const gsl::not_null<Variables<DgPackageFieldTags>*> lower_packaged_data,
71  const gsl::not_null<Variables<DgPackageFieldTags>*> upper_packaged_data,
72  const size_t logical_dimension_to_operate_in, const Element<Dim>& element,
73  const Mesh<Dim>& subcell_volume_mesh,
74  const std::unordered_map<
77  boost::hash<std::pair<Direction<Dim>, ElementId<Dim>>>>&
78  mortar_data) noexcept {
79  const Direction<Dim> upper_direction{logical_dimension_to_operate_in,
80  Side::Upper};
81  const Direction<Dim> lower_direction{logical_dimension_to_operate_in,
82  Side::Lower};
83  const bool has_upper_neighbor = element.neighbors().contains(upper_direction);
84  const bool has_lower_neighbor = element.neighbors().contains(lower_direction);
85  const std::pair upper_mortar_id =
86  has_upper_neighbor
87  ? std::pair{upper_direction,
88  *element.neighbors().at(upper_direction).begin()}
90  const std::pair lower_mortar_id =
91  has_lower_neighbor
92  ? std::pair{lower_direction,
93  *element.neighbors().at(lower_direction).begin()}
95 
96  Index<Dim> subcell_extents_with_faces = subcell_volume_mesh.extents();
97  ++subcell_extents_with_faces[logical_dimension_to_operate_in];
98  const Mesh<Dim - 1>& subcell_face_mesh =
99  subcell_volume_mesh.slice_away(logical_dimension_to_operate_in);
100 
101  const auto project_dg_data_to_subcells =
102  [logical_dimension_to_operate_in, &subcell_extents_with_faces,
103  &subcell_face_mesh](const gsl::not_null<Variables<DgPackageFieldTags>*>
104  subcell_packaged_data,
105  const size_t subcell_index,
106  const Mesh<Dim - 1>& neighbor_face_mesh,
107  const std::vector<double>& neighbor_data) noexcept {
108  const double* slice_data = neighbor_data.data();
109  // Warning: projected_data can't be inside the `if constexpr` since that
110  // would lead to a dangling pointer.
111  std::vector<double> projected_data{};
112  if constexpr (Dim > 1) {
113  projected_data.resize(
114  Variables<DgPackageFieldTags>::number_of_independent_components *
115  subcell_face_mesh.number_of_grid_points());
116  evolution::dg::subcell::fd::detail::project_impl(
117  gsl::make_span(projected_data.data(), projected_data.size()),
118  gsl::make_span(neighbor_data.data(), neighbor_data.size()),
119  neighbor_face_mesh, subcell_face_mesh.extents());
120  slice_data = projected_data.data();
121  } else {
122  (void)subcell_face_mesh;
123  (void)projected_data;
124  }
125  const size_t volume_grid_points = subcell_extents_with_faces.product();
126  const size_t slice_grid_points =
127  subcell_extents_with_faces
128  .slice_away(logical_dimension_to_operate_in)
129  .product();
130  double* const volume_data = subcell_packaged_data->data();
131  for (SliceIterator si(subcell_extents_with_faces,
132  logical_dimension_to_operate_in, subcell_index);
133  si; ++si) {
134  for (size_t i = 0;
135  i <
136  Variables<DgPackageFieldTags>::number_of_independent_components;
137  ++i) {
138  // clang-tidy: do not use pointer arithmetic
139  volume_data[si.volume_offset() +
140  i * volume_grid_points] = // NOLINT
141  slice_data[si.slice_offset() +
142  i * slice_grid_points]; // NOLINT
143  }
144  }
145  };
146 
147  // Project DG data to the subcells
148  if (auto neighbor_mortar_data_it = mortar_data.find(upper_mortar_id);
149  has_upper_neighbor and neighbor_mortar_data_it != mortar_data.end()) {
150  if (neighbor_mortar_data_it->second.neighbor_mortar_data().has_value()) {
151  project_dg_data_to_subcells(
152  upper_packaged_data,
153  subcell_extents_with_faces[logical_dimension_to_operate_in] - 1,
154  neighbor_mortar_data_it->second.neighbor_mortar_data()->first,
155  neighbor_mortar_data_it->second.neighbor_mortar_data()->second);
156  }
157  if constexpr (OverwriteInternalMortarData) {
158  if (neighbor_mortar_data_it->second.local_mortar_data().has_value()) {
159  project_dg_data_to_subcells(
160  lower_packaged_data,
161  subcell_extents_with_faces[logical_dimension_to_operate_in] - 1,
162  neighbor_mortar_data_it->second.local_mortar_data()->first,
163  neighbor_mortar_data_it->second.local_mortar_data()->second);
164  }
165  }
166  }
167  if (auto neighbor_mortar_data_it = mortar_data.find(lower_mortar_id);
168  has_lower_neighbor and neighbor_mortar_data_it != mortar_data.end()) {
169  if (neighbor_mortar_data_it->second.neighbor_mortar_data().has_value()) {
170  project_dg_data_to_subcells(
171  lower_packaged_data, 0,
172  neighbor_mortar_data_it->second.neighbor_mortar_data()->first,
173  neighbor_mortar_data_it->second.neighbor_mortar_data()->second);
174  }
175  if constexpr (OverwriteInternalMortarData) {
176  if (neighbor_mortar_data_it->second.local_mortar_data().has_value()) {
177  project_dg_data_to_subcells(
178  upper_packaged_data, 0,
179  neighbor_mortar_data_it->second.local_mortar_data()->first,
180  neighbor_mortar_data_it->second.local_mortar_data()->second);
181  }
182  }
183  }
184 }
185 } // namespace evolution::dg::subcell
utility
std::pair
evolution::dg::subcell::correct_package_data
void correct_package_data(const gsl::not_null< Variables< DgPackageFieldTags > * > lower_packaged_data, const gsl::not_null< Variables< DgPackageFieldTags > * > upper_packaged_data, const size_t logical_dimension_to_operate_in, const Element< Dim > &element, const Mesh< Dim > &subcell_volume_mesh, const std::unordered_map< std::pair< Direction< Dim >, ElementId< Dim >>, evolution::dg::MortarData< Dim >, boost::hash< std::pair< Direction< Dim >, ElementId< Dim >>>> &mortar_data) noexcept
Project the DG package data to the subcells. Data received from a neighboring element doing DG is alw...
Definition: CorrectPackagedData.hpp:69
std::vector< double >
Index
Definition: Index.hpp:31
algorithm
Direction< Dim >
Element
Definition: Element.hpp:29
ElementId< Dim >
ElementId.hpp
evolution::dg::MortarData
Data on the mortar used to compute the boundary correction for the DG scheme.
Definition: MortarData.hpp:50
cstddef
Assert.hpp
SliceIterator
Iterate over a (dim-1)-dimensional slice.
Definition: SliceIterator.hpp:23
Index.hpp
Variables.hpp
Element.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
gsl::make_span
constexpr span< ElementType > make_span(ElementType *ptr, typename span< ElementType >::index_type count)
Definition: Gsl.hpp:801
evolution::dg::subcell
Implementation of a generic finite volume/conservative finite difference subcell limiter.
Definition: Actions.hpp:6
Gsl.hpp
Tensor.hpp
evolution::dg::subcell::slice_data
DirectionMap< Dim, std::vector< double > > slice_data(const Variables< TagList > &volume_subcell_vars, const Index< Dim > &subcell_extents, const size_t number_of_ghost_points, const DirectionMap< Dim, bool > &directions_to_slice) noexcept
Slice the subcell variables needed for neighbors to perform reconstruction.
Definition: SliceData.hpp:45
Direction.hpp
Mesh::slice_away
Mesh< Dim - 1 > slice_away(size_t d) const noexcept
Returns a Mesh with dimension d removed (zero-indexed).
unordered_map
TMPL.hpp
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13