ElementCenteredSubdomainData.hpp
1 // 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 
14 #include "ParallelAlgorithms/LinearSolver/Schwarz/OverlapHelpers.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>
36  static constexpr size_t volume_dim = Dim;
37  using ElementData = Variables<TagsList>;
38  using OverlapData = ElementData;
40  using const_iterator =
42 
43  ElementCenteredSubdomainData() = default;
46  const ElementCenteredSubdomainData&) noexcept = default;
48  default;
50  ElementCenteredSubdomainData&&) noexcept = default;
51  ~ElementCenteredSubdomainData() noexcept = default;
52 
54  const size_t element_num_points) noexcept
55  : element_data{element_num_points} {}
56 
57  template <typename UsedForSizeTagsList>
58  void destructive_resize(
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 
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  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  iterator begin() noexcept { return {this}; }
90  iterator end() noexcept { return {}; }
91  const_iterator begin() const noexcept { return {this}; }
92  const_iterator end() const noexcept { return {}; }
93  const_iterator cbegin() const noexcept { return begin(); }
94  const_iterator cend() const noexcept { return end(); }
95 
96  void pup(PUP::er& p) noexcept { // NOLINT
97  p | element_data;
98  p | overlap_data;
99  }
100 
101  template <typename RhsTagsList>
102  ElementCenteredSubdomainData& operator+=(
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  ElementCenteredSubdomainData& operator-=(
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  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  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  ElementData element_data{};
142  OverlapMap<Dim, OverlapData> overlap_data{};
143 };
144 
145 template <size_t Dim, typename LhsTagsList, typename RhsTagsList>
146 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 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 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 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 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 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 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 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_id(d);
233  const auto rhs_segment_id = rhs.second.segment_id(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>
265  private:
266  using PtrType =
267  tmpl::conditional_t<Const,
270 
271  public:
272  using difference_type = ptrdiff_t;
273  using value_type = double;
274  using pointer = value_type*;
275  using reference = value_type&;
277 
278  /// Construct begin state
279  ElementCenteredSubdomainDataIterator(PtrType data) noexcept : data_(data) {
280  overlap_ids_ = detail::ordered_overlap_ids(data_->overlap_data);
281  reset();
282  }
283 
284  void reset() noexcept {
285  overlap_index_ = (data_->element_data.size() == 0 and overlap_ids_.empty())
287  : 0;
288  data_index_ = 0;
289  }
290 
291  /// Construct end state
293  overlap_index_ = std::numeric_limits<size_t>::max();
294  data_index_ = 0;
295  }
296 
297  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  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  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  friend bool operator!=(
331  const ElementCenteredSubdomainDataIterator& lhs,
332  const ElementCenteredSubdomainDataIterator& rhs) noexcept {
333  return not(lhs == rhs);
334  }
335 
336  PtrType data_;
337  std::vector<OverlapId<Dim>> overlap_ids_;
338  size_t overlap_index_;
339  size_t data_index_;
340 };
341 
342 } // namespace LinearSolver::Schwarz
343 
345 
346 template <size_t Dim, typename LhsTagsList, typename RhsTagsList>
348  Schwarz::ElementCenteredSubdomainData<Dim, LhsTagsList>,
349  Schwarz::ElementCenteredSubdomainData<Dim, RhsTagsList>> {
350  static double apply(
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>
368  LinearSolver::Schwarz::ElementCenteredSubdomainData<Dim, TagsListOut>,
370  using SubdomainDataIn =
372  using SubdomainDataOut =
375  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
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
InnerProduct.hpp
std::rel_ops::operator!=
T operator!=(T... args)
std::vector
MakeWithValueImpls::MakeWithValueImpl
Definition: MakeWithValue.hpp:64
LinearSolver::Schwarz::ElementCenteredSubdomainData
Data on an element-centered subdomain.
Definition: ElementCenteredSubdomainData.hpp:35
std::forward_iterator_tag
tuple
MakeWithValueImpls::MakeWithValueImpl::apply
static R apply(const T &input, const ValueType value) noexcept
The default implementation uses number_of_points and MakeWithSize.
Definition: MakeWithValue.hpp:67
LinearSolver::Schwarz::ElementCenteredSubdomainDataIterator::ElementCenteredSubdomainDataIterator
ElementCenteredSubdomainDataIterator(PtrType data) noexcept
Construct begin state.
Definition: ElementCenteredSubdomainData.hpp:279
LinearSolver::Schwarz
Items related to the Schwarz linear solver.
Definition: CommunicateOverlapFields.hpp:37
LinearSolver::Schwarz::ElementCenteredSubdomainDataIterator
Iterate over LinearSolver::Schwarz::ElementCenteredSubdomainData
Definition: ElementCenteredSubdomainData.hpp:264
LinearSolver::inner_product
double inner_product(const Lhs &lhs, const Rhs &rhs) noexcept
The local part of the Euclidean inner product on the vector space w.r.t. which the addition and scala...
Definition: InnerProduct.hpp:71
SPECTRE_ALWAYS_INLINE
#define SPECTRE_ALWAYS_INLINE
Definition: ForceInline.hpp:16
std::ostream
cstddef
MakeWithValueImpls
Definition: DenseVector.hpp:61
MakeWithValue.hpp
operator*
auto operator*(const TensorExpression< T1, typename T1::type, typename T1::symmetry, typename T1::index_list, ArgsList1 > &t1, const TensorExpression< T2, typename T2::type, typename T2::symmetry, typename T2::index_list, ArgsList2 > &t2)
Returns the tensor expression representing the product of two tensor expressions.
Definition: Product.hpp:223
cpp2b::accumulate
constexpr T accumulate(InputIt first, InputIt last, T init)
Definition: Numeric.hpp:33
LinearSolver::InnerProductImpls
Definition: InnerProduct.hpp:20
Variables.hpp
LinearSolver
Functionality for solving linear systems of equations.
Definition: ExplicitInverse.hpp:20
std::numeric_limits::max
T max(T... args)
LinearSolver::InnerProductImpls::InnerProductImpl
The inner product between any types that have a dot product.
Definition: InnerProduct.hpp:24
FixedHashMap
A hash table with a compile-time specified maximum size and ability to efficiently handle perfect has...
Definition: FixedHashMap.hpp:82
ostream
std::numeric_limits
std::data
T data(T... args)
LinearSolver::Schwarz::ElementCenteredSubdomainDataIterator::ElementCenteredSubdomainDataIterator
ElementCenteredSubdomainDataIterator() noexcept
Construct end state.
Definition: ElementCenteredSubdomainData.hpp:292
TMPL.hpp
string