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 <ostream>
8 : #include <pup.h>
9 : #include <string>
10 : #include <tuple>
11 :
12 : #include "DataStructures/Variables.hpp"
13 : #include "NumericalAlgorithms/LinearSolver/InnerProduct.hpp"
14 : #include "ParallelAlgorithms/LinearSolver/Schwarz/OverlapHelpers.hpp"
15 : #include "Utilities/MakeWithValue.hpp"
16 : #include "Utilities/TMPL.hpp"
17 :
18 : namespace LinearSolver::Schwarz {
19 :
20 : /// \cond
21 : template <bool Const, size_t Dim, typename TagsList>
22 : struct ElementCenteredSubdomainDataIterator;
23 : /// \endcond
24 :
25 : /*!
26 : * \brief Data on an element-centered subdomain
27 : *
28 : * An element-centered subdomain consists of a central element and overlap
29 : * regions with all neighboring elements. This class holds data on such a
30 : * subdomain. It supports vector space operations (addition and scalar
31 : * multiplication) and an inner product, which allows the use of this data type
32 : * with linear solvers (see e.g. `LinearSolver::Serial::Gmres`).
33 : */
34 : template <size_t Dim, typename TagsList>
35 1 : struct ElementCenteredSubdomainData {
36 0 : static constexpr size_t volume_dim = Dim;
37 0 : using ElementData = Variables<TagsList>;
38 0 : using OverlapData = ElementData;
39 0 : using iterator = ElementCenteredSubdomainDataIterator<false, Dim, TagsList>;
40 0 : using const_iterator =
41 : ElementCenteredSubdomainDataIterator<true, Dim, TagsList>;
42 :
43 0 : ElementCenteredSubdomainData() = default;
44 0 : ElementCenteredSubdomainData(const ElementCenteredSubdomainData&) = default;
45 0 : ElementCenteredSubdomainData& operator=(
46 : const ElementCenteredSubdomainData&) noexcept = default;
47 0 : ElementCenteredSubdomainData(ElementCenteredSubdomainData&&) noexcept =
48 : default;
49 0 : ElementCenteredSubdomainData& operator=(
50 : ElementCenteredSubdomainData&&) noexcept = default;
51 0 : ~ElementCenteredSubdomainData() noexcept = default;
52 :
53 0 : explicit ElementCenteredSubdomainData(
54 : const size_t element_num_points) noexcept
55 : : element_data{element_num_points} {}
56 :
57 : template <typename UsedForSizeTagsList>
58 0 : void destructive_resize(
59 : const ElementCenteredSubdomainData<Dim, UsedForSizeTagsList>&
60 : used_for_size) noexcept {
61 : if (UNLIKELY(element_data.number_of_grid_points() !=
62 : used_for_size.element_data.number_of_grid_points())) {
63 : element_data.initialize(
64 : used_for_size.element_data.number_of_grid_points());
65 : }
66 : for (const auto& [overlap_id, used_for_overlap_size] :
67 : used_for_size.overlap_data) {
68 : if (UNLIKELY(overlap_data[overlap_id].number_of_grid_points() !=
69 : used_for_overlap_size.number_of_grid_points())) {
70 : overlap_data[overlap_id].initialize(
71 : used_for_overlap_size.number_of_grid_points());
72 : }
73 : }
74 : }
75 :
76 0 : ElementCenteredSubdomainData(
77 : Variables<TagsList> local_element_data,
78 : OverlapMap<Dim, Variables<TagsList>> local_overlap_data) noexcept
79 : : element_data{std::move(local_element_data)},
80 : overlap_data{std::move(local_overlap_data)} {}
81 :
82 0 : size_t size() const noexcept {
83 : return std::accumulate(
84 : overlap_data.begin(), overlap_data.end(), element_data.size(),
85 : [](const size_t size, const auto& overlap_id_and_data) noexcept {
86 : return size + overlap_id_and_data.second.size();
87 : });
88 : }
89 0 : iterator begin() noexcept { return {this}; }
90 0 : iterator end() noexcept { return {}; }
91 0 : const_iterator begin() const noexcept { return {this}; }
92 0 : const_iterator end() const noexcept { return {}; }
93 0 : const_iterator cbegin() const noexcept { return begin(); }
94 0 : const_iterator cend() const noexcept { return end(); }
95 :
96 0 : void pup(PUP::er& p) noexcept { // NOLINT
97 : p | element_data;
98 : p | overlap_data;
99 : }
100 :
101 : template <typename RhsTagsList>
102 0 : ElementCenteredSubdomainData& operator+=(
103 : const ElementCenteredSubdomainData<Dim, RhsTagsList>& rhs) noexcept {
104 : element_data += rhs.element_data;
105 : for (auto& [overlap_id, data] : overlap_data) {
106 : data += rhs.overlap_data.at(overlap_id);
107 : }
108 : return *this;
109 : }
110 :
111 : template <typename RhsTagsList>
112 0 : ElementCenteredSubdomainData& operator-=(
113 : const ElementCenteredSubdomainData<Dim, RhsTagsList>& rhs) noexcept {
114 : element_data -= rhs.element_data;
115 : for (auto& [overlap_id, data] : overlap_data) {
116 : data -= rhs.overlap_data.at(overlap_id);
117 : }
118 : return *this;
119 : }
120 :
121 0 : ElementCenteredSubdomainData& operator*=(const double scalar) noexcept {
122 : element_data *= scalar;
123 : for (auto& [overlap_id, data] : overlap_data) {
124 : data *= scalar;
125 : // Silence unused-variable warning on GCC 7
126 : (void)overlap_id;
127 : }
128 : return *this;
129 : }
130 :
131 0 : ElementCenteredSubdomainData& operator/=(const double scalar) noexcept {
132 : element_data /= scalar;
133 : for (auto& [overlap_id, data] : overlap_data) {
134 : data /= scalar;
135 : // Silence unused-variable warning on GCC 7
136 : (void)overlap_id;
137 : }
138 : return *this;
139 : }
140 :
141 0 : ElementData element_data{};
142 0 : OverlapMap<Dim, OverlapData> overlap_data{};
143 : };
144 :
145 : template <size_t Dim, typename LhsTagsList, typename RhsTagsList>
146 0 : decltype(auto) operator+(
147 : ElementCenteredSubdomainData<Dim, LhsTagsList> lhs,
148 : const ElementCenteredSubdomainData<Dim, RhsTagsList>& rhs) noexcept {
149 : lhs += rhs;
150 : return lhs;
151 : }
152 :
153 : template <size_t Dim, typename LhsTagsList, typename RhsTagsList>
154 0 : decltype(auto) operator-(
155 : ElementCenteredSubdomainData<Dim, LhsTagsList> lhs,
156 : const ElementCenteredSubdomainData<Dim, RhsTagsList>& rhs) noexcept {
157 : lhs -= rhs;
158 : return lhs;
159 : }
160 :
161 : template <size_t Dim, typename TagsList>
162 0 : decltype(auto) operator*(
163 : const double scalar,
164 : ElementCenteredSubdomainData<Dim, TagsList> data) noexcept {
165 : data *= scalar;
166 : return data;
167 : }
168 :
169 : template <size_t Dim, typename TagsList>
170 0 : decltype(auto) operator*(ElementCenteredSubdomainData<Dim, TagsList> data,
171 : const double scalar) noexcept {
172 : data *= scalar;
173 : return data;
174 : }
175 :
176 : template <size_t Dim, typename TagsList>
177 0 : decltype(auto) operator/(ElementCenteredSubdomainData<Dim, TagsList> data,
178 : const double scalar) noexcept {
179 : data /= scalar;
180 : return data;
181 : }
182 :
183 : template <size_t Dim, typename TagsList>
184 0 : std::ostream& operator<<(std::ostream& os,
185 : const ElementCenteredSubdomainData<Dim, TagsList>&
186 : subdomain_data) noexcept {
187 : os << "Element data:\n"
188 : << subdomain_data.element_data << "\nOverlap data:\n"
189 : << subdomain_data.overlap_data;
190 : return os;
191 : }
192 :
193 : template <size_t Dim, typename TagsList>
194 0 : bool operator==(
195 : const ElementCenteredSubdomainData<Dim, TagsList>& lhs,
196 : const ElementCenteredSubdomainData<Dim, TagsList>& rhs) noexcept {
197 : return lhs.element_data == rhs.element_data and
198 : lhs.overlap_data == rhs.overlap_data;
199 : }
200 :
201 : template <size_t Dim, typename TagsList>
202 0 : bool operator!=(
203 : const ElementCenteredSubdomainData<Dim, TagsList>& lhs,
204 : const ElementCenteredSubdomainData<Dim, TagsList>& rhs) noexcept {
205 : return not(lhs == rhs);
206 : }
207 :
208 : namespace detail {
209 : // Defines a consistent ordering of overlap IDs
210 : template <size_t Dim, typename ValueType>
211 : std::vector<OverlapId<Dim>> ordered_overlap_ids(
212 : const OverlapMap<Dim, ValueType>& overlap_map) noexcept {
213 : std::vector<OverlapId<Dim>> overlap_ids{};
214 : overlap_ids.reserve(overlap_map.size());
215 : std::transform(overlap_map.begin(), overlap_map.end(),
216 : std::back_inserter(overlap_ids),
217 : [](const auto& overlap_id_and_value) noexcept {
218 : return overlap_id_and_value.first;
219 : });
220 : std::sort(overlap_ids.begin(), overlap_ids.end(),
221 : [](const OverlapId<Dim>& lhs, const OverlapId<Dim>& rhs) noexcept {
222 : if (lhs.first.axis() != rhs.first.axis()) {
223 : return lhs.first.axis() < rhs.first.axis();
224 : }
225 : if (lhs.first.side() != rhs.first.side()) {
226 : return lhs.first.side() < rhs.first.side();
227 : }
228 : if (lhs.second.block_id() != rhs.second.block_id()) {
229 : return lhs.second.block_id() < rhs.second.block_id();
230 : }
231 : for (size_t d = 0; d < Dim; ++d) {
232 : const auto& lhs_segment_id = lhs.second.segment_ids().at(d);
233 : const auto& rhs_segment_id = rhs.second.segment_ids().at(d);
234 : if (lhs_segment_id.refinement_level() !=
235 : rhs_segment_id.refinement_level()) {
236 : return lhs_segment_id.refinement_level() <
237 : rhs_segment_id.refinement_level();
238 : }
239 : if (lhs_segment_id.index() != rhs_segment_id.index()) {
240 : return lhs_segment_id.index() < rhs_segment_id.index();
241 : }
242 : }
243 : return false;
244 : });
245 : return overlap_ids;
246 : }
247 : } // namespace detail
248 :
249 : /*!
250 : * \brief Iterate over `LinearSolver::Schwarz::ElementCenteredSubdomainData`
251 : *
252 : * This iterator guarantees that it steps through the data in the same order as
253 : * long as these conditions are satisfied:
254 : *
255 : * - The set of overlap IDs in the `overlap_data` doesn't change
256 : * - The extents of the `element_data` and the `overlap_data doesn't change
257 : *
258 : * Iterating requires sorting the overlap IDs. If you find this impacts
259 : * performance, be advised to implement the internal data storage in
260 : * `ElementCenteredSubdomainData` so it stores its data contiguously, e.g. by
261 : * implementing non-owning variables.
262 : */
263 : template <bool Const, size_t Dim, typename TagsList>
264 1 : struct ElementCenteredSubdomainDataIterator {
265 : private:
266 0 : using PtrType =
267 : tmpl::conditional_t<Const,
268 : const ElementCenteredSubdomainData<Dim, TagsList>*,
269 : ElementCenteredSubdomainData<Dim, TagsList>*>;
270 :
271 : public:
272 0 : using difference_type = ptrdiff_t;
273 0 : using value_type = double;
274 0 : using pointer = value_type*;
275 0 : using reference = value_type&;
276 0 : using iterator_category = std::forward_iterator_tag;
277 :
278 : /// Construct begin state
279 1 : ElementCenteredSubdomainDataIterator(PtrType data) noexcept : data_(data) {
280 : overlap_ids_ = detail::ordered_overlap_ids(data_->overlap_data);
281 : reset();
282 : }
283 :
284 0 : void reset() noexcept {
285 : overlap_index_ = (data_->element_data.size() == 0 and overlap_ids_.empty())
286 : ? std::numeric_limits<size_t>::max()
287 : : 0;
288 : data_index_ = 0;
289 : }
290 :
291 : /// Construct end state
292 1 : ElementCenteredSubdomainDataIterator() noexcept {
293 : overlap_index_ = std::numeric_limits<size_t>::max();
294 : data_index_ = 0;
295 : }
296 :
297 0 : ElementCenteredSubdomainDataIterator& operator++() noexcept {
298 : ++data_index_;
299 : if (data_index_ ==
300 : (overlap_index_ == 0
301 : ? data_->element_data
302 : : data_->overlap_data.at(overlap_ids_[overlap_index_ - 1]))
303 : .size()) {
304 : ++overlap_index_;
305 : data_index_ = 0;
306 : }
307 : if (overlap_index_ == overlap_ids_.size() + 1) {
308 : overlap_index_ = std::numeric_limits<size_t>::max();
309 : }
310 : return *this;
311 : }
312 :
313 0 : tmpl::conditional_t<Const, double, double&> operator*() const noexcept {
314 : if (overlap_index_ == 0) {
315 : return data_->element_data.data()[data_index_];
316 : } else {
317 : return data_->overlap_data.at(overlap_ids_[overlap_index_ - 1])
318 : .data()[data_index_];
319 : }
320 : }
321 :
322 : private:
323 0 : friend bool operator==(
324 : const ElementCenteredSubdomainDataIterator& lhs,
325 : const ElementCenteredSubdomainDataIterator& rhs) noexcept {
326 : return lhs.overlap_index_ == rhs.overlap_index_ and
327 : lhs.data_index_ == rhs.data_index_;
328 : }
329 :
330 0 : friend bool operator!=(
331 : const ElementCenteredSubdomainDataIterator& lhs,
332 : const ElementCenteredSubdomainDataIterator& rhs) noexcept {
333 : return not(lhs == rhs);
334 : }
335 :
336 0 : PtrType data_;
337 0 : std::vector<OverlapId<Dim>> overlap_ids_;
338 0 : size_t overlap_index_;
339 0 : size_t data_index_;
340 : };
341 :
342 : } // namespace LinearSolver::Schwarz
343 :
344 : namespace LinearSolver::InnerProductImpls {
345 :
346 : template <size_t Dim, typename LhsTagsList, typename RhsTagsList>
347 : struct InnerProductImpl<
348 : Schwarz::ElementCenteredSubdomainData<Dim, LhsTagsList>,
349 0 : Schwarz::ElementCenteredSubdomainData<Dim, RhsTagsList>> {
350 0 : static double apply(
351 : const Schwarz::ElementCenteredSubdomainData<Dim, LhsTagsList>& lhs,
352 : const Schwarz::ElementCenteredSubdomainData<Dim, RhsTagsList>&
353 : rhs) noexcept {
354 : double result = inner_product(lhs.element_data, rhs.element_data);
355 : for (const auto& [overlap_id, lhs_data] : lhs.overlap_data) {
356 : result += inner_product(lhs_data, rhs.overlap_data.at(overlap_id));
357 : }
358 : return result;
359 : }
360 : };
361 :
362 : } // namespace LinearSolver::InnerProductImpls
363 :
364 : namespace MakeWithValueImpls {
365 :
366 : template <size_t Dim, typename TagsListOut, typename TagsListIn>
367 : struct MakeWithValueImpl<
368 : LinearSolver::Schwarz::ElementCenteredSubdomainData<Dim, TagsListOut>,
369 0 : LinearSolver::Schwarz::ElementCenteredSubdomainData<Dim, TagsListIn>> {
370 0 : using SubdomainDataIn =
371 : LinearSolver::Schwarz::ElementCenteredSubdomainData<Dim, TagsListIn>;
372 0 : using SubdomainDataOut =
373 : LinearSolver::Schwarz::ElementCenteredSubdomainData<Dim, TagsListOut>;
374 : static SPECTRE_ALWAYS_INLINE SubdomainDataOut
375 0 : apply(const SubdomainDataIn& input, const double value) noexcept {
376 : SubdomainDataOut output{};
377 : output.element_data =
378 : make_with_value<typename SubdomainDataOut::ElementData>(
379 : input.element_data, value);
380 : for (const auto& [overlap_id, input_data] : input.overlap_data) {
381 : output.overlap_data.emplace(
382 : overlap_id, make_with_value<typename SubdomainDataOut::OverlapData>(
383 : input_data, value));
384 : }
385 : return output;
386 : }
387 : };
388 :
389 : } // namespace MakeWithValueImpls
|
