OverlapHelpers.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <boost/functional/hash.hpp>
7 #include <cstddef>
8 #include <limits>
9 #include <tuple>
10 
11 #include "DataStructures/DataVector.hpp"
12 #include "DataStructures/FixedHashMap.hpp"
13 #include "DataStructures/Index.hpp"
18 #include "Domain/Structure/MaxNumberOfNeighbors.hpp"
19 #include "ErrorHandling/Assert.hpp"
20 #include "Utilities/ContainerHelpers.hpp"
21 
22 namespace LinearSolver::Schwarz {
23 
24 /// Identifies a subdomain region that overlaps with another element
25 template <size_t Dim>
27 
28 /// Data structure that can store the `ValueType` on each possible overlap of an
29 /// element-centered subdomain with its neighbors. Overlaps are identified by
30 /// their `OverlapId`.
31 template <size_t Dim, typename ValueType>
32 using OverlapMap =
34  boost::hash<OverlapId<Dim>>>;
35 
36 /*!
37  * \brief The number of points that an overlap extends into the `volume_extent`
38  *
39  * In a dimension where an element has `volume_extent` points, the overlap
40  * extent is the largest number under these constraints:
41  *
42  * - It is at most `max_overlap`.
43  * - It is smaller than the `volume_extent`.
44  *
45  * This means the overlap extent is always smaller than the `volume_extent`. The
46  * reason for this constraint is that we define the _width_ of the overlap as
47  * the element-logical coordinate distance from the face of the element to the
48  * first collocation point _outside_ the overlap extent. Therefore, even an
49  * overlap region that covers the full element in width does not include the
50  * collocation point on the opposite side of the element.
51  *
52  * Here's a few notes on the definition of the overlap extent and width:
53  *
54  * - A typical smooth weighting function goes to zero at the overlap width, so
55  * if the grid points located at the overlap width were included in the
56  * subdomain, their solutions would not contribute to the weighted sum of
57  * subdomain solutions.
58  * - Defining the overlap width as the distance to the first point _outside_ the
59  * overlap extent makes it non-zero even for a single point of overlap into a
60  * Gauss-Lobatto grid (which has points located at the element face).
61  * - Boundary contributions for many (but not all) discontinuous Galerkin
62  * schemes on Gauss-Lobatto grids are limited to the grid points on the
63  * element face, e.g. for a DG operator that is pre-multiplied by the mass
64  * matrix, or one where boundary contributions are lifted using the diagonal
65  * mass-matrix approximation. Not including the grid points facing away from
66  * the subdomain in the overlap allows to ignore that face altogether in the
67  * subdomain operator.
68  */
69 size_t overlap_extent(size_t volume_extent, size_t max_overlap) noexcept;
70 
71 /*!
72  * \brief Total number of grid points in an overlap region that extends
73  * `overlap_extent` points into the `volume_extents` from either side in the
74  * `overlap_dimension`
75  *
76  * The overlap region has `overlap_extent` points in the `overlap_dimension`,
77  * and `volume_extents` points in the other dimensions. The number of grid
78  * points returned by this function is the product of these extents.
79  */
80 template <size_t Dim>
81 size_t overlap_num_points(const Index<Dim>& volume_extents,
82  size_t overlap_extent,
83  size_t overlap_dimension) noexcept;
84 
85 /*!
86  * \brief Width of an overlap extending `overlap_extent` points into the
87  * `collocation_points` from either side.
88  *
89  * The "width" of an overlap is the element-logical coordinate distance from the
90  * element boundary to the first collocation point outside the overlap region in
91  * the overlap dimension, i.e. the dimension perpendicular to the element face.
92  * See `LinearSolver::Schwarz::overlap_extent` for details.
93  *
94  * This function assumes the `collocation_points` are mirrored around 0.
95  */
96 double overlap_width(size_t overlap_extent,
97  const DataVector& collocation_points) noexcept;
98 
99 /*!
100  * \brief Iterate over grid points in a region that extends partially into the
101  * volume
102  *
103  * Here's an example how to use this iterator:
104  *
105  * \snippet Test_OverlapHelpers.cpp overlap_iterator
106  */
108  public:
109  template <size_t Dim>
110  OverlapIterator(const Index<Dim>& volume_extents, size_t overlap_extent,
111  const Direction<Dim>& direction) noexcept;
112 
113  explicit operator bool() const noexcept;
114 
115  OverlapIterator& operator++();
116 
117  /// Offset into a DataVector that holds full volume data
118  size_t volume_offset() const noexcept;
119 
120  /// Offset into a DataVector that holds data only on the overlap region
121  size_t overlap_offset() const noexcept;
122 
123  void reset() noexcept;
124 
125  private:
126  size_t size_ = std::numeric_limits<size_t>::max();
127  size_t num_slices_ = std::numeric_limits<size_t>::max();
128  size_t stride_ = std::numeric_limits<size_t>::max();
129  size_t stride_count_ = std::numeric_limits<size_t>::max();
130  size_t jump_ = std::numeric_limits<size_t>::max();
131  size_t initial_offset_ = std::numeric_limits<size_t>::max();
132  size_t volume_offset_ = std::numeric_limits<size_t>::max();
133  size_t overlap_offset_ = std::numeric_limits<size_t>::max();
134 };
135 
136 // @{
137 /// The part of the tensor data that lies within the overlap region
138 template <size_t Dim, typename DataType, typename... TensorStructure>
139 void data_on_overlap(const gsl::not_null<Tensor<DataType, TensorStructure...>*>
140  restricted_tensor,
141  const Tensor<DataType, TensorStructure...>& tensor,
142  const Index<Dim>& volume_extents,
143  const size_t overlap_extent,
144  const Direction<Dim>& direction) noexcept {
145  for (OverlapIterator overlap_iterator{volume_extents, overlap_extent,
146  direction};
147  overlap_iterator; ++overlap_iterator) {
148  for (size_t tensor_component = 0; tensor_component < tensor.size();
149  ++tensor_component) {
150  (*restricted_tensor)[tensor_component][overlap_iterator
151  .overlap_offset()] =
152  tensor[tensor_component][overlap_iterator.volume_offset()];
153  }
154  }
155 }
156 
157 template <size_t Dim, typename DataType, typename... TensorStructure>
158 Tensor<DataType, TensorStructure...> data_on_overlap(
159  const Tensor<DataType, TensorStructure...>& tensor,
160  const Index<Dim>& volume_extents, const size_t overlap_extent,
161  const Direction<Dim>& direction) noexcept {
162  Tensor<DataType, TensorStructure...> restricted_tensor{overlap_num_points(
163  volume_extents, overlap_extent, direction.dimension())};
164  data_on_overlap(make_not_null(&restricted_tensor), tensor, volume_extents,
165  overlap_extent, direction);
166  return restricted_tensor;
167 }
168 
169 namespace detail {
170 template <size_t Dim>
171 void data_on_overlap_impl(double* overlap_data, const double* volume_data,
172  size_t num_components,
173  const Index<Dim>& volume_extents,
174  size_t overlap_extent,
175  const Direction<Dim>& direction) noexcept;
176 } // namespace detail
177 
178 template <size_t Dim, typename OverlapTagsList, typename VolumeTagsList>
179 void data_on_overlap(
180  const gsl::not_null<Variables<OverlapTagsList>*> overlap_data,
181  const Variables<VolumeTagsList>& volume_data,
182  const Index<Dim>& volume_extents, const size_t overlap_extent,
183  const Direction<Dim>& direction) noexcept {
184  constexpr size_t num_components =
185  Variables<VolumeTagsList>::number_of_independent_components;
186  ASSERT(volume_data.number_of_grid_points() == volume_extents.product(),
187  "volume_data has wrong number of grid points. Expected "
188  << volume_extents.product() << ", got "
189  << volume_data.number_of_grid_points());
190  ASSERT(overlap_data->number_of_grid_points() ==
191  overlap_num_points(volume_extents, overlap_extent,
192  direction.dimension()),
193  "overlap_data has wrong number of grid points. Expected "
194  << overlap_num_points(volume_extents, overlap_extent,
195  direction.dimension())
196  << ", got " << overlap_data->number_of_grid_points());
197  detail::data_on_overlap_impl(overlap_data->data(), volume_data.data(),
198  num_components, volume_extents, overlap_extent,
199  direction);
200 }
201 
202 template <size_t Dim, typename TagsList>
203 Variables<TagsList> data_on_overlap(const Variables<TagsList>& volume_data,
204  const Index<Dim>& volume_extents,
205  const size_t overlap_extent,
206  const Direction<Dim>& direction) noexcept {
207  Variables<TagsList> overlap_data{overlap_num_points(
208  volume_extents, overlap_extent, direction.dimension())};
209  data_on_overlap(make_not_null(&overlap_data), volume_data, volume_extents,
210  overlap_extent, direction);
211  return overlap_data;
212 }
213 // @}
214 
215 namespace detail {
216 template <size_t Dim>
217 void add_overlap_data_impl(double* volume_data, const double* overlap_data,
218  size_t num_components,
219  const Index<Dim>& volume_extents,
220  size_t overlap_extent,
221  const Direction<Dim>& direction) noexcept;
222 } // namespace detail
223 
224 /// Add the `overlap_data` to the `volume_data`
225 template <size_t Dim, typename VolumeTagsList, typename OverlapTagsList>
227  const gsl::not_null<Variables<VolumeTagsList>*> volume_data,
228  const Variables<OverlapTagsList>& overlap_data,
229  const Index<Dim>& volume_extents, const size_t overlap_extent,
230  const Direction<Dim>& direction) noexcept {
231  constexpr size_t num_components =
232  Variables<VolumeTagsList>::number_of_independent_components;
233  ASSERT(volume_data->number_of_grid_points() == volume_extents.product(),
234  "volume_data has wrong number of grid points. Expected "
235  << volume_extents.product() << ", got "
236  << volume_data->number_of_grid_points());
237  ASSERT(overlap_data.number_of_grid_points() ==
238  overlap_num_points(volume_extents, overlap_extent,
239  direction.dimension()),
240  "overlap_data has wrong number of grid points. Expected "
241  << overlap_num_points(volume_extents, overlap_extent,
242  direction.dimension())
243  << ", got " << overlap_data.number_of_grid_points());
244  detail::add_overlap_data_impl(volume_data->data(), overlap_data.data(),
245  num_components, volume_extents, overlap_extent,
246  direction);
247 }
248 
249 // @{
250 /// Extend the overlap data to the full mesh by filling it with zeros outside
251 /// the overlap region
252 template <size_t Dim, typename ExtendedTagsList, typename OverlapTagsList>
254  const gsl::not_null<Variables<ExtendedTagsList>*> extended_data,
255  const Variables<OverlapTagsList>& overlap_data,
256  const Index<Dim>& volume_extents, const size_t overlap_extent,
257  const Direction<Dim>& direction) noexcept {
258  *extended_data = Variables<ExtendedTagsList>{volume_extents.product(), 0.};
259  add_overlap_data(extended_data, overlap_data, volume_extents, overlap_extent,
260  direction);
261 }
262 
263 template <size_t Dim, typename TagsList>
264 Variables<TagsList> extended_overlap_data(
265  const Variables<TagsList>& overlap_data, const Index<Dim>& volume_extents,
266  const size_t overlap_extent, const Direction<Dim>& direction) noexcept {
267  Variables<TagsList> extended_data{volume_extents.product()};
268  extended_overlap_data(make_not_null(&extended_data), overlap_data,
269  volume_extents, overlap_extent, direction);
270  return extended_data;
271 }
272 // @}
273 
274 } // namespace LinearSolver::Schwarz
LinearSolver::Schwarz::OverlapIterator::volume_offset
size_t volume_offset() const noexcept
Offset into a DataVector that holds full volume data.
Definition: OverlapHelpers.cpp:104
maximum_number_of_neighbors
constexpr size_t maximum_number_of_neighbors(const size_t dim)
Definition: MaxNumberOfNeighbors.hpp:15
LinearSolver::Schwarz::overlap_extent
size_t overlap_extent(const size_t volume_extent, const size_t max_overlap) noexcept
The number of points that an overlap extends into the volume_extent
Definition: OverlapHelpers.cpp:21
std::pair
tuple
Index
Definition: Index.hpp:31
LinearSolver::Schwarz::OverlapIterator::overlap_offset
size_t overlap_offset() const noexcept
Offset into a DataVector that holds data only on the overlap region.
Definition: OverlapHelpers.cpp:108
Direction
Definition: Direction.hpp:23
LinearSolver::Schwarz
Items related to the Schwarz linear solver.
Definition: CommunicateOverlapFields.hpp:36
ElementId< Dim >
ElementId.hpp
LinearSolver::Schwarz::OverlapIterator
Iterate over grid points in a region that extends partially into the volume.
Definition: OverlapHelpers.hpp:107
LinearSolver::Schwarz::add_overlap_data
void add_overlap_data(const gsl::not_null< Variables< VolumeTagsList > * > volume_data, const Variables< OverlapTagsList > &overlap_data, const Index< Dim > &volume_extents, const size_t overlap_extent, const Direction< Dim > &direction) noexcept
Add the overlap_data to the volume_data
Definition: OverlapHelpers.hpp:226
LinearSolver::Schwarz::overlap_width
double overlap_width(const size_t overlap_extent, const DataVector &collocation_points) noexcept
Width of an overlap extending overlap_extent points into the collocation_points from either side.
Definition: OverlapHelpers.cpp:49
cstddef
Assert.hpp
LinearSolver::Schwarz::data_on_overlap
void data_on_overlap(const gsl::not_null< Tensor< DataType, TensorStructure... > * > restricted_tensor, const Tensor< DataType, TensorStructure... > &tensor, const Index< Dim > &volume_extents, const size_t overlap_extent, const Direction< Dim > &direction) noexcept
The part of the tensor data that lies within the overlap region.
Definition: OverlapHelpers.hpp:139
Index.hpp
LinearSolver::Schwarz::overlap_num_points
size_t overlap_num_points(const Index< Dim > &volume_extents, const size_t overlap_extent, const size_t overlap_dimension) noexcept
Total number of grid points in an overlap region that extends overlap_extent points into the volume_e...
Definition: OverlapHelpers.cpp:41
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:42
LinearSolver::Schwarz::extended_overlap_data
void extended_overlap_data(const gsl::not_null< Variables< ExtendedTagsList > * > extended_data, const Variables< OverlapTagsList > &overlap_data, const Index< Dim > &volume_extents, const size_t overlap_extent, const Direction< Dim > &direction) noexcept
Extend the overlap data to the full mesh by filling it with zeros outside the overlap region.
Definition: OverlapHelpers.hpp:253
Spectral::collocation_points
const DataVector & collocation_points(const Mesh< 1 > &mesh) noexcept
Collocation points for a one-dimensional mesh.
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:51
Variables.hpp
limits
Tensor.hpp
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
std::numeric_limits::max
T max(T... args)
FixedHashMap
A hash table with a compile-time specified maximum size and ability to efficiently handle perfect has...
Definition: FixedHashMap.hpp:81
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183