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 
17 namespace LinearSolver::Schwarz {
18 
19 /*!
20  * \brief Data on an element-centered subdomain
21  *
22  * An element-centered subdomain consists of a central element and overlap
23  * regions with all neighboring elements. This class holds data on such a
24  * subdomain. It supports vector space operations (addition and scalar
25  * multiplication) and an inner product, which allows the use of this data type
26  * with linear solvers (see e.g. `LinearSolver::Serial::Gmres`).
27  */
28 template <size_t Dim, typename TagsList>
30  static constexpr size_t volume_dim = Dim;
31  using ElementData = Variables<TagsList>;
32  using OverlapData = ElementData;
33 
34  ElementCenteredSubdomainData() = default;
37  const ElementCenteredSubdomainData&) noexcept = default;
39  default;
41  ElementCenteredSubdomainData&&) noexcept = default;
42  ~ElementCenteredSubdomainData() noexcept = default;
43 
45  const size_t element_num_points) noexcept
46  : element_data{element_num_points} {}
47 
48  template <typename UsedForSizeTagsList>
49  void destructive_resize(
51  used_for_size) noexcept {
52  if (UNLIKELY(element_data.number_of_grid_points() !=
53  used_for_size.element_data.number_of_grid_points())) {
54  element_data.initialize(
55  used_for_size.element_data.number_of_grid_points());
56  }
57  for (const auto& [overlap_id, used_for_overlap_size] :
58  used_for_size.overlap_data) {
59  if (UNLIKELY(overlap_data[overlap_id].number_of_grid_points() !=
60  used_for_overlap_size.number_of_grid_points())) {
61  overlap_data[overlap_id].initialize(
62  used_for_overlap_size.number_of_grid_points());
63  }
64  }
65  }
66 
68  Variables<TagsList> local_element_data,
69  OverlapMap<Dim, Variables<TagsList>> local_overlap_data) noexcept
70  : element_data{std::move(local_element_data)},
71  overlap_data{std::move(local_overlap_data)} {}
72 
73  void pup(PUP::er& p) noexcept { // NOLINT
74  p | element_data;
75  p | overlap_data;
76  }
77 
78  template <typename RhsTagsList>
79  ElementCenteredSubdomainData& operator+=(
81  element_data += rhs.element_data;
82  for (auto& [overlap_id, data] : overlap_data) {
83  data += rhs.overlap_data.at(overlap_id);
84  }
85  return *this;
86  }
87 
88  template <typename RhsTagsList>
89  ElementCenteredSubdomainData& operator-=(
91  element_data -= rhs.element_data;
92  for (auto& [overlap_id, data] : overlap_data) {
93  data -= rhs.overlap_data.at(overlap_id);
94  }
95  return *this;
96  }
97 
98  ElementCenteredSubdomainData& operator*=(const double scalar) noexcept {
99  element_data *= scalar;
100  for (auto& [overlap_id, data] : overlap_data) {
101  data *= scalar;
102  // Silence unused-variable warning on GCC 7
103  (void)overlap_id;
104  }
105  return *this;
106  }
107 
108  ElementCenteredSubdomainData& operator/=(const double scalar) noexcept {
109  element_data /= scalar;
110  for (auto& [overlap_id, data] : overlap_data) {
111  data /= scalar;
112  // Silence unused-variable warning on GCC 7
113  (void)overlap_id;
114  }
115  return *this;
116  }
117 
118  ElementData element_data{};
119  OverlapMap<Dim, OverlapData> overlap_data{};
120 };
121 
122 template <size_t Dim, typename LhsTagsList, typename RhsTagsList>
123 decltype(auto) operator+(
124  ElementCenteredSubdomainData<Dim, LhsTagsList> lhs,
125  const ElementCenteredSubdomainData<Dim, RhsTagsList>& rhs) noexcept {
126  lhs += rhs;
127  return lhs;
128 }
129 
130 template <size_t Dim, typename LhsTagsList, typename RhsTagsList>
131 decltype(auto) operator-(
132  ElementCenteredSubdomainData<Dim, LhsTagsList> lhs,
133  const ElementCenteredSubdomainData<Dim, RhsTagsList>& rhs) noexcept {
134  lhs -= rhs;
135  return lhs;
136 }
137 
138 template <size_t Dim, typename TagsList>
139 decltype(auto) operator*(
140  const double scalar,
141  ElementCenteredSubdomainData<Dim, TagsList> data) noexcept {
142  data *= scalar;
143  return data;
144 }
145 
146 template <size_t Dim, typename TagsList>
147 decltype(auto) operator*(ElementCenteredSubdomainData<Dim, TagsList> data,
148  const double scalar) noexcept {
149  data *= scalar;
150  return data;
151 }
152 
153 template <size_t Dim, typename TagsList>
154 decltype(auto) operator/(ElementCenteredSubdomainData<Dim, TagsList> data,
155  const double scalar) noexcept {
156  data /= scalar;
157  return data;
158 }
159 
160 template <size_t Dim, typename TagsList>
161 std::ostream& operator<<(std::ostream& os,
162  const ElementCenteredSubdomainData<Dim, TagsList>&
163  subdomain_data) noexcept {
164  os << "Element data:\n"
165  << subdomain_data.element_data << "\nOverlap data:\n"
166  << subdomain_data.overlap_data;
167  return os;
168 }
169 
170 template <size_t Dim, typename TagsList>
171 bool operator==(
172  const ElementCenteredSubdomainData<Dim, TagsList>& lhs,
173  const ElementCenteredSubdomainData<Dim, TagsList>& rhs) noexcept {
174  return lhs.element_data == rhs.element_data and
175  lhs.overlap_data == rhs.overlap_data;
176 }
177 
178 template <size_t Dim, typename TagsList>
179 bool operator!=(
180  const ElementCenteredSubdomainData<Dim, TagsList>& lhs,
181  const ElementCenteredSubdomainData<Dim, TagsList>& rhs) noexcept {
182  return not(lhs == rhs);
183 }
184 
185 } // namespace LinearSolver::Schwarz
186 
188 
189 template <size_t Dim, typename LhsTagsList, typename RhsTagsList>
191  Schwarz::ElementCenteredSubdomainData<Dim, LhsTagsList>,
192  Schwarz::ElementCenteredSubdomainData<Dim, RhsTagsList>> {
193  static double apply(
196  rhs) noexcept {
197  double result = inner_product(lhs.element_data, rhs.element_data);
198  for (const auto& [overlap_id, lhs_data] : lhs.overlap_data) {
199  result += inner_product(lhs_data, rhs.overlap_data.at(overlap_id));
200  }
201  return result;
202  }
203 };
204 
205 } // namespace LinearSolver::InnerProductImpls
206 
207 namespace MakeWithValueImpls {
208 
209 template <size_t Dim, typename TagsListOut, typename TagsListIn>
211  LinearSolver::Schwarz::ElementCenteredSubdomainData<Dim, TagsListOut>,
213  using SubdomainDataIn =
215  using SubdomainDataOut =
218  apply(const SubdomainDataIn& input, const double value) noexcept {
219  SubdomainDataOut output{};
220  output.element_data =
221  make_with_value<typename SubdomainDataOut::ElementData>(
222  input.element_data, value);
223  for (const auto& [overlap_id, input_data] : input.overlap_data) {
224  output.overlap_data.emplace(
225  overlap_id, make_with_value<typename SubdomainDataOut::OverlapData>(
226  input_data, value));
227  }
228  return output;
229  }
230 };
231 
232 } // namespace MakeWithValueImpls
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
InnerProduct.hpp
std::rel_ops::operator!=
T operator!=(T... args)
MakeWithValueImpls::MakeWithValueImpl
Definition: MakeWithValue.hpp:64
LinearSolver::Schwarz::ElementCenteredSubdomainData
Data on an element-centered subdomain.
Definition: ElementCenteredSubdomainData.hpp:29
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
Items related to the Schwarz linear solver.
Definition: CommunicateOverlapFields.hpp:36
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
LinearSolver::InnerProductImpls
Definition: InnerProduct.hpp:20
Variables.hpp
LinearSolver
Functionality for solving linear systems of equations.
Definition: Gmres.cpp:16
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::data
T data(T... args)
string