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 double 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 double 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 0 : decltype(auto) operator*(const double scalar,
214 : ElementCenteredSubdomainData<Dim, TagsList> data) {
215 : data *= scalar;
216 : return data;
217 : }
218 :
219 : template <size_t Dim, typename TagsList>
220 0 : decltype(auto) operator*(ElementCenteredSubdomainData<Dim, TagsList> data,
221 : const double scalar) {
222 : data *= scalar;
223 : return data;
224 : }
225 :
226 : template <size_t Dim, typename TagsList>
227 0 : decltype(auto) operator/(ElementCenteredSubdomainData<Dim, TagsList> data,
228 : const double scalar) {
229 : data /= scalar;
230 : return data;
231 : }
232 :
233 : template <size_t Dim, typename TagsList>
234 0 : std::ostream& operator<<(
235 : std::ostream& os,
236 : const ElementCenteredSubdomainData<Dim, TagsList>& subdomain_data) {
237 : os << "Element data:\n"
238 : << subdomain_data.element_data << "\nOverlap data:\n"
239 : << subdomain_data.overlap_data;
240 : return os;
241 : }
242 :
243 : template <size_t Dim, typename TagsList>
244 0 : bool operator==(const ElementCenteredSubdomainData<Dim, TagsList>& lhs,
245 : const ElementCenteredSubdomainData<Dim, TagsList>& rhs) {
246 : return lhs.element_data == rhs.element_data and
247 : lhs.overlap_data == rhs.overlap_data;
248 : }
249 :
250 : template <size_t Dim, typename TagsList>
251 0 : bool operator!=(const ElementCenteredSubdomainData<Dim, TagsList>& lhs,
252 : const ElementCenteredSubdomainData<Dim, TagsList>& rhs) {
253 : return not(lhs == rhs);
254 : }
255 :
256 : /*!
257 : * \brief Iterate over `LinearSolver::Schwarz::ElementCenteredSubdomainData`
258 : *
259 : * This iterator guarantees that it steps through the data in the same order as
260 : * long as these conditions are satisfied:
261 : *
262 : * - The set of overlap IDs in the `overlap_data` doesn't change
263 : * - The extents of the `element_data` and the `overlap_data doesn't change
264 : *
265 : * Iterating requires sorting the overlap IDs. If you find this impacts
266 : * performance, be advised to implement the internal data storage in
267 : * `ElementCenteredSubdomainData` so it stores its data contiguously, e.g. by
268 : * implementing non-owning variables.
269 : */
270 : template <bool Const, size_t Dim, typename TagsList>
271 1 : struct ElementCenteredSubdomainDataIterator {
272 : private:
273 0 : using PtrType =
274 : tmpl::conditional_t<Const,
275 : const ElementCenteredSubdomainData<Dim, TagsList>*,
276 : ElementCenteredSubdomainData<Dim, TagsList>*>;
277 0 : using VarsPtrType = tmpl::conditional_t<Const, const Variables<TagsList>*,
278 : Variables<TagsList>*>;
279 :
280 : public:
281 0 : using difference_type = ptrdiff_t;
282 0 : using value_type = double;
283 0 : using pointer = value_type*;
284 0 : using reference = value_type&;
285 0 : using iterator_category = std::forward_iterator_tag;
286 :
287 : /// Construct begin state
288 1 : ElementCenteredSubdomainDataIterator(PtrType data) : data_(data) {
289 : reset();
290 : }
291 :
292 0 : void reset() {
293 : overlap_index_ = (data_->element_data.size() == 0 and
294 : data_->ordered_overlap_ids_.empty())
295 : ? std::numeric_limits<size_t>::max()
296 : : 0;
297 : update_vars_ref();
298 : data_index_ = 0;
299 : }
300 :
301 : /// Construct end state
302 1 : ElementCenteredSubdomainDataIterator() {
303 : overlap_index_ = std::numeric_limits<size_t>::max();
304 : data_index_ = 0;
305 : }
306 :
307 0 : ElementCenteredSubdomainDataIterator& operator++() {
308 : ++data_index_;
309 : if (data_index_ == vars_ref_->size()) {
310 : ++overlap_index_;
311 : update_vars_ref();
312 : data_index_ = 0;
313 : }
314 : if (overlap_index_ == data_->ordered_overlap_ids_.size() + 1) {
315 : overlap_index_ = std::numeric_limits<size_t>::max();
316 : }
317 : return *this;
318 : }
319 :
320 0 : tmpl::conditional_t<Const, double, double&> operator*() const {
321 : return vars_ref_->data()[data_index_];
322 : }
323 :
324 : private:
325 0 : void update_vars_ref() {
326 : // The random-access operation is relatively slow because it computes a
327 : // hash, so it's important for performance to avoid repeating it at every
328 : // data index. Instead, we update this reference only when the overlap index
329 : // changes.
330 : if (overlap_index_ > data_->ordered_overlap_ids_.size()) {
331 : return;
332 : }
333 : vars_ref_ = overlap_index_ == 0
334 : ? &data_->element_data
335 : : &data_->overlap_data.at(
336 : data_->ordered_overlap_ids_[overlap_index_ - 1]);
337 : }
338 :
339 0 : friend bool operator==(const ElementCenteredSubdomainDataIterator& lhs,
340 : const ElementCenteredSubdomainDataIterator& rhs) {
341 : return lhs.overlap_index_ == rhs.overlap_index_ and
342 : lhs.data_index_ == rhs.data_index_;
343 : }
344 :
345 0 : friend bool operator!=(const ElementCenteredSubdomainDataIterator& lhs,
346 : const ElementCenteredSubdomainDataIterator& rhs) {
347 : return not(lhs == rhs);
348 : }
349 :
350 0 : PtrType data_;
351 0 : size_t overlap_index_;
352 0 : size_t data_index_;
353 0 : VarsPtrType vars_ref_ = nullptr;
354 : };
355 :
356 : } // namespace LinearSolver::Schwarz
357 :
358 : namespace LinearSolver::InnerProductImpls {
359 :
360 : template <size_t Dim, typename LhsTagsList, typename RhsTagsList>
361 0 : struct InnerProductImpl<
362 : Schwarz::ElementCenteredSubdomainData<Dim, LhsTagsList>,
363 : Schwarz::ElementCenteredSubdomainData<Dim, RhsTagsList>> {
364 0 : static double apply(
365 : const Schwarz::ElementCenteredSubdomainData<Dim, LhsTagsList>& lhs,
366 : const Schwarz::ElementCenteredSubdomainData<Dim, RhsTagsList>& rhs) {
367 : double result = inner_product(lhs.element_data, rhs.element_data);
368 : for (const auto& [overlap_id, lhs_data] : lhs.overlap_data) {
369 : result += inner_product(lhs_data, rhs.overlap_data.at(overlap_id));
370 : }
371 : return result;
372 : }
373 : };
374 :
375 : } // namespace LinearSolver::InnerProductImpls
376 :
377 : namespace MakeWithValueImpls {
378 :
379 : template <size_t Dim, typename TagsListOut, typename TagsListIn>
380 0 : struct MakeWithValueImpl<
381 : LinearSolver::Schwarz::ElementCenteredSubdomainData<Dim, TagsListOut>,
382 : LinearSolver::Schwarz::ElementCenteredSubdomainData<Dim, TagsListIn>> {
383 0 : using SubdomainDataIn =
384 : LinearSolver::Schwarz::ElementCenteredSubdomainData<Dim, TagsListIn>;
385 0 : using SubdomainDataOut =
386 : LinearSolver::Schwarz::ElementCenteredSubdomainData<Dim, TagsListOut>;
387 : static SPECTRE_ALWAYS_INLINE SubdomainDataOut
388 0 : apply(const SubdomainDataIn& input, const double value) {
389 : SubdomainDataOut output{};
390 : output.element_data =
391 : make_with_value<typename SubdomainDataOut::ElementData>(
392 : input.element_data, value);
393 : for (const auto& [overlap_id, input_data] : input.overlap_data) {
394 : output.overlap_data.emplace(
395 : overlap_id, make_with_value<typename SubdomainDataOut::OverlapData>(
396 : input_data, value));
397 : }
398 : return output;
399 : }
400 : };
401 :
402 : } // namespace MakeWithValueImpls
|
