GridFormat 0.4.0
I/O-Library for grid-like data structures
Loading...
Searching...
No Matches
parallel.hpp
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2022-2023 Dennis Gläser <dennis.glaeser@iws.uni-stuttgart.de>
2// SPDX-License-Identifier: MIT
8#ifndef GRIDFORMAT_VTK_PARALLEL_HPP_
9#define GRIDFORMAT_VTK_PARALLEL_HPP_
10
11#include <array>
12#include <vector>
13#include <string>
14#include <utility>
15#include <concepts>
16#include <algorithm>
17#include <iterator>
18#include <numeric>
19#include <tuple>
20#include <cmath>
21
22#include <gridformat/common/math.hpp>
23#include <gridformat/common/exceptions.hpp>
24#include <gridformat/common/concepts.hpp>
25#include <gridformat/common/ranges.hpp>
27#include <gridformat/grid/grid.hpp>
28
29#include <gridformat/parallel/communication.hpp>
30#include <gridformat/parallel/helpers.hpp>
31
32#include <gridformat/xml/element.hpp>
34
35namespace GridFormat::PVTK {
36
38 std::string piece_basefilename(const std::string& par_filename, int rank) {
39 const std::string base_name = par_filename.substr(0, par_filename.find_last_of("."));
40 return base_name + "-" + std::to_string(rank);
41}
42
44template<typename Encoder, typename DataFormat>
46 public:
47 PDataArrayHelper(const Encoder& e,
48 const DataFormat& df,
49 XMLElement& element)
50 : _encoder(e)
51 , _data_format(df)
52 , _element(element)
53 {}
54
55 void add(const std::string& name, const Field& field) {
56 // vtk always uses 3d, this assumes that the field
57 // is transformed accordingly in the piece writers
58 static constexpr std::size_t vtk_space_dim = 3;
59 const auto layout = field.layout();
60 std::size_t ncomps = std::pow(vtk_space_dim, field.layout().dimension() - 1);
61
62 // leave higher-dimensional vectors untouched
63 if (layout.dimension() == 2 && layout.extent(1) > vtk_space_dim)
64 ncomps = layout.extent(1);
65
66 XMLElement& arr = _element.add_child("PDataArray");
67 arr.set_attribute("Name", name);
68 arr.set_attribute("type", VTK::attribute_name(field.precision()));
69 arr.set_attribute("format", VTK::data_format_name(_encoder, _data_format));
70 arr.set_attribute("NumberOfComponents", ncomps);
71 }
72
73 private:
74 const Encoder& _encoder;
75 const DataFormat& _data_format;
76 XMLElement& _element;
77};
78
79template<std::size_t dim>
81 public:
82 using Location = std::array<std::size_t, dim>;
83
84 explicit StructuredGridMapper(std::vector<Location>&& map)
85 : _map{std::move(map)}
86 {}
87
88 const auto& location(std::integral auto rank) const {
89 return _map[rank];
90 }
91
92 std::ranges::view auto ranks_below(const Location& loc, unsigned direction) const {
93 return std::views::filter(
94 std::views::iota(std::size_t{0}, _map.size()),
95 [&, loc=loc, direction=direction] (const std::size_t rank) {
96 const auto& rank_loc = _map[rank];
97 for (unsigned dir = 0; dir < dim; ++dir)
98 if (dir != direction && rank_loc[dir] != loc[dir]) {
99 return false;
100 }
101 return rank_loc[direction] < loc[direction];
102 }
103 );
104 }
105
106 private:
107 std::vector<Location> _map;
108};
109
112template<typename T, std::size_t dim>
114 public:
115 using Origin = std::array<T, dim>;
116 using Basis = std::array<Origin, dim>;
117
118 explicit StructuredGridMapperHelper(std::integral auto ranks,
119 const Basis& basis,
120 T default_epsilon = 1e-6)
121 : _basis(basis)
122 , _origins(ranks)
123 , _set(ranks, false)
124 , _default_epsilon{default_epsilon} {
125 std::ranges::fill(_reverse, false);
126 }
127
128 void reverse(int direction) {
129 _reverse[direction] = !_reverse[direction];
130 }
131
132 void set_origin_for(int rank, Origin origin) {
133 if (_set[rank])
134 throw ValueError("Origin for given rank already set");
135 _origins[rank] = std::move(origin);
136 _set[rank] = true;
137 }
138
139 auto make_mapper() const {
140 if (_origins.empty())
141 throw InvalidState("No origins have been set");
142
143 std::vector<std::array<std::size_t, dim>> map(_origins.size());
144 for (unsigned dir = 0; dir < dim; ++dir) {
145 auto ordinates = _get_ordinates(dir);
146 const auto eps = _epsilon(ordinates);
147 _sort_ordinates(ordinates, eps, _reverse[dir]);
148
149 if (ordinates.size() > 1) {
150 for (unsigned rank = 0; rank < _origins.size(); ++rank) {
151 const auto rank_origin_ordinate = _get_ordinate(_origins[rank], dir);
152 auto it = std::ranges::find_if(ordinates, [&] (const T o) {
153 using std::abs;
154 return abs(o - rank_origin_ordinate) < eps;
155 });
156 if (it == ordinates.end())
157 throw InvalidState("Could not determine rank ordinate");
158 map[rank][dir] = std::distance(ordinates.begin(), it);
159 }
160 } else {
161 for (unsigned rank = 0; rank < _origins.size(); ++rank)
162 map[rank][dir] = 0;
163 }
164 }
165 return StructuredGridMapper<dim>{std::move(map)};
166 }
167
168 Origin compute_origin() const {
169 if (_origins.empty())
170 throw InvalidState("No origins have been set");
171
172 Origin result;
173 for (unsigned dir = 0; dir < dim; ++dir) {
174 auto ordinates = _get_ordinates(dir);
175 _sort_ordinates(ordinates, _reverse[dir]);
176 result[dir] = ordinates[0];
177 }
178 return result;
179 }
180
181 private:
182 auto _get_ordinates(unsigned int axis) const {
183 std::vector<T> result;
184 result.reserve(_origins.size());
185 std::ranges::copy(
186 std::views::transform(_origins, [&] (const Origin& o) {
187 return _get_ordinate(o, axis);
188 }),
189 std::back_inserter(result)
190 );
191 return result;
192 }
193
194 auto _get_ordinate(const Origin& o, unsigned int axis) const {
195 return dot_product(o, _basis[axis]);
196 }
197
198 void _sort_ordinates(std::vector<T>& ordinates, bool reverse) const {
199 if (ordinates.size() > 1)
200 _sort_ordinates(ordinates, _epsilon(ordinates), reverse);
201 }
202
203 void _sort_ordinates(std::vector<T>& ordinates, const T& eps, bool reverse) const {
204 if (ordinates.size() > 1) {
205 auto [it, _] = Ranges::sort_and_unique(ordinates, {}, [&] (T a, T b) {
206 using std::abs;
207 return abs(a - b) < eps;
208 });
209 ordinates.erase(it, ordinates.end());
210 if (reverse)
211 std::ranges::reverse(ordinates);
212 }
213 }
214
215 T _epsilon(const std::vector<T>& ordinates) const {
216 const auto max = *std::ranges::max_element(ordinates);
217 const auto min = *std::ranges::min_element(ordinates);
218 const auto size = max - min;
219
220 if (size > _default_epsilon) {
221 std::vector<T> dx; dx.reserve(ordinates.size());
222 std::adjacent_difference(ordinates.begin(), ordinates.end(), std::back_inserter(dx));
223 std::ranges::for_each(dx, [] (auto& v) { v = std::abs(v); });
224 std::ranges::sort(dx);
225 for (T _dx : dx)
226 if (_dx > 1e-8*size)
227 return _dx*0.1;
228 }
229
230 return _default_epsilon;
231 }
232
233 const Basis& _basis;
234 std::vector<Origin> _origins;
235 std::vector<bool> _set;
236 std::array<bool, dim> _reverse;
237 T _default_epsilon;
238};
239
240
241template<typename Communicator>
243 public:
244 StructuredParallelGridHelper(const Communicator& comm, int root_rank = 0)
245 : _comm{comm}
246 , _root_rank{root_rank}
247 {}
248
249 template<Concepts::Scalar CT,
250 std::size_t dim,
251 Concepts::StaticallySizedMDRange<2> B = std::array<std::array<CT, dim>, dim>>
252 auto compute_extents_and_origin(const std::vector<CT>& all_origins,
253 const std::vector<std::size_t>& all_extents,
254 const std::array<bool, dim>& is_negative_axis,
255 const B& basis = GridDetail::standard_basis<CT, dim>()) const {
256 const auto num_ranks = Parallel::size(_comm);
257 std::vector<std::array<std::size_t, dim>> pieces_begin(num_ranks);
258 std::vector<std::array<std::size_t, dim>> pieces_end(num_ranks);
259 std::array<std::size_t, dim> whole_extent;
260 std::array<CT, dim> origin;
261
262 if (Parallel::rank(_comm) == _root_rank) {
263 const auto default_epsilon = 1e-6*std::ranges::max(all_origins | std::views::transform([] (CT v) {
264 using std::abs;
265 return abs(v);
266 }));
267 const auto mapper_helper = _make_mapper_helper(basis, all_origins, is_negative_axis, default_epsilon);
268 const auto rank_mapper = mapper_helper.make_mapper();
269 origin = mapper_helper.compute_origin();
270
271 for (unsigned dir = 0; dir < dim; ++dir) {
272 for (int rank = 0; rank < num_ranks; ++rank) {
273 auto ranks_below = rank_mapper.ranks_below(rank_mapper.location(rank), dir);
274 const std::size_t offset = std::accumulate(
275 std::ranges::begin(ranks_below),
276 std::ranges::end(ranks_below),
277 std::size_t{0},
278 [&] (const std::size_t current, int r) {
279 return current + Parallel::access_gathered<dim>(all_extents, _comm, {dir, r});
280 }
281 );
282 pieces_begin[rank][dir] = offset;
283 pieces_end[rank][dir] = offset + Parallel::access_gathered<dim>(all_extents, _comm, {dir, rank});
284 }
285
286 whole_extent[dir] = (*std::max_element(
287 pieces_end.begin(), pieces_end.end(),
288 [&] (const auto& a1, const auto& a2) { return a1[dir] < a2[dir]; }
289 ))[dir];
290 }
291 }
292
293 return std::make_tuple(
294 std::move(pieces_begin),
295 std::move(pieces_end),
296 std::move(whole_extent),
297 std::move(origin)
298 );
299 }
300
301 private:
302 template<Concepts::StaticallySizedMDRange<2> Basis, Concepts::Scalar CT, std::size_t dim>
303 auto _make_mapper_helper(const Basis& basis,
304 const std::vector<CT>& all_origins,
305 const std::array<bool, dim>& is_negative_axis,
306 CT default_eps) const {
307 StructuredGridMapperHelper<CT, dim> helper(Parallel::size(_comm), basis, default_eps);
308 std::ranges::for_each(Parallel::ranks(_comm), [&] (int rank) {
309 helper.set_origin_for(rank, Parallel::access_gathered<dim>(all_origins, _comm, rank));
310 });
311
312 for (unsigned dir = 0; dir < dim; ++dir)
313 if (is_negative_axis[dir])
314 helper.reverse(dir);
315
316 return helper;
317 }
318
319 const Communicator& _comm;
320 int _root_rank;
321};
322
323} // namespace GridFormat::PVTK
324
325#endif // GRIDFORMAT_VTK_PARALLEL_HPP_
Helper functions to get the VTK-specific names of things.
Abstract interface for fields of values that is used by writers/readers to store fields.
Definition: field.hpp:37
DynamicPrecision precision() const
Return the precision of the scalar field values.
Definition: field.hpp:56
MDLayout layout() const
Return the layout of this field.
Definition: field.hpp:51
Helper to add a PDataArray child to an xml element.
Definition: parallel.hpp:45
Definition: parallel.hpp:80
std::string piece_basefilename(const std::string &par_filename, int rank)
Return the piece filename (w/o extension) for the given rank.
Definition: parallel.hpp:38