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