Line data Source code
1 0 : // 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"
14 : #include "DataStructures/Tensor/Tensor.hpp"
15 : #include "DataStructures/Variables.hpp"
16 : #include "Domain/Structure/Direction.hpp"
17 : #include "Domain/Structure/ElementId.hpp"
18 : #include "Domain/Structure/MaxNumberOfNeighbors.hpp"
19 : #include "Utilities/ContainerHelpers.hpp"
20 : #include "Utilities/ErrorHandling/Assert.hpp"
21 :
22 : namespace LinearSolver::Schwarz {
23 :
24 : /// Identifies a subdomain region that overlaps with another element
25 : template <size_t Dim>
26 1 : using OverlapId = std::pair<Direction<Dim>, ElementId<Dim>>;
27 :
28 : /// Data structure that can store the `ValueType` on each possible overlap of an
29 : /// elementcentered subdomain with its neighbors. Overlaps are identified by
30 : /// their `OverlapId`.
31 : template <size_t Dim, typename ValueType>
32 1 : using OverlapMap =
33 : FixedHashMap<maximum_number_of_neighbors(Dim), OverlapId<Dim>, ValueType,
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 elementlogical 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 nonzero even for a single point of overlap into a
60 : * GaussLobatto grid (which has points located at the element face).
61 : *  Boundary contributions for many (but not all) discontinuous Galerkin
62 : * schemes on GaussLobatto grids are limited to the grid points on the
63 : * element face, e.g. for a DG operator that is premultiplied by the mass
64 : * matrix, or one where boundary contributions are lifted using the diagonal
65 : * massmatrix 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 1 : 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 1 : 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 elementlogical 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 1 : 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 : */
107 1 : class OverlapIterator {
108 : public:
109 : template <size_t Dim>
110 0 : OverlapIterator(const Index<Dim>& volume_extents, size_t overlap_extent,
111 : const Direction<Dim>& direction) noexcept;
112 :
113 0 : explicit operator bool() const noexcept;
114 :
115 0 : OverlapIterator& operator++();
116 :
117 : /// Offset into a DataVector that holds full volume data
118 1 : size_t volume_offset() const noexcept;
119 :
120 : /// Offset into a DataVector that holds data only on the overlap region
121 1 : size_t overlap_offset() const noexcept;
122 :
123 0 : void reset() noexcept;
124 :
125 : private:
126 0 : size_t size_ = std::numeric_limits<size_t>::max();
127 0 : size_t num_slices_ = std::numeric_limits<size_t>::max();
128 0 : size_t stride_ = std::numeric_limits<size_t>::max();
129 0 : size_t stride_count_ = std::numeric_limits<size_t>::max();
130 0 : size_t jump_ = std::numeric_limits<size_t>::max();
131 0 : size_t initial_offset_ = std::numeric_limits<size_t>::max();
132 0 : size_t volume_offset_ = std::numeric_limits<size_t>::max();
133 0 : 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 1 : 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 0 : 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 0 : 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 0 : 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>
226 1 : void add_overlap_data(
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>
253 1 : void extended_overlap_data(
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 0 : 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
