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