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