GCC Code Coverage Report

Directory: gridformat/
File: gridformat/vtk/parallel.hpp
Date: 2024-11-20 14:41:59
Exec Total Coverage
Lines: 149 155 96.1%
Functions: 51 71 71.8%
Branches: 131 222 59.0%
Line Branch Exec Source
1 // SPDX-FileCopyrightText: 2022-2023 Dennis Gläser <dennis.glaeser@iws.uni-stuttgart.de>
2 // SPDX-License-Identifier: MIT
3 /*!
4 * \file
5 * \ingroup VTK
6 * \brief Helper function for writing parallel VTK files
7 */
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>
26 #include <gridformat/common/field.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>
33 #include <gridformat/vtk/attributes.hpp>
34
35 namespace GridFormat::PVTK {
36
37 //! Return the piece filename (w/o extension) for the given rank
38 12740 std::string piece_basefilename(const std::string& par_filename, int rank) {
39
1/2
✓ Branch 2 taken 12740 times.
✗ Branch 3 not taken.
 12740 const std::string base_name = par_filename.substr(0, par_filename.find_last_of("."));
40
2/4
✓ Branch 2 taken 12740 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12740 times.
✗ Branch 6 not taken.
 25480 return base_name + "-" + std::to_string(rank);
41 12740 }
42
43 //! Helper to add a PDataArray child to an xml element
44 template<typename Encoder, typename DataFormat>
45 class PDataArrayHelper {
46 public:
47 8480 PDataArrayHelper(const Encoder& e,
48 const DataFormat& df,
49 XMLElement& element)
50 8480 : _encoder(e)
51 8480 , _data_format(df)
52 8480 , _element(element)
53 8480 {}
54
55 46968 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
1/2
✓ Branch 1 taken 23484 times.
✗ Branch 2 not taken.
 46968 const auto layout = field.layout();
60
2/4
✓ Branch 1 taken 23484 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23484 times.
✗ Branch 5 not taken.
 46968 std::size_t ncomps = std::pow(vtk_space_dim, field.layout().dimension() - 1);
61
62 // leave higher-dimensional vectors untouched
63
8/10
✓ Branch 1 taken 23484 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 7910 times.
✓ Branch 4 taken 15574 times.
✓ Branch 6 taken 7910 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 144 times.
✓ Branch 9 taken 7766 times.
✓ Branch 10 taken 144 times.
✓ Branch 11 taken 23340 times.
 46968 if (layout.dimension() == 2 && layout.extent(1) > vtk_space_dim)
64
1/2
✓ Branch 1 taken 144 times.
✗ Branch 2 not taken.
 288 ncomps = layout.extent(1);
65
66
2/4
✓ Branch 1 taken 23484 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23484 times.
✗ Branch 5 not taken.
 93936 XMLElement& arr = _element.add_child("PDataArray");
67
2/4
✓ Branch 1 taken 23484 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23484 times.
✗ Branch 5 not taken.
 93936 arr.set_attribute("Name", name);
68
4/8
✓ Branch 1 taken 23484 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23484 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 23484 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 23484 times.
✗ Branch 11 not taken.
 140904 arr.set_attribute("type", VTK::attribute_name(field.precision()));
69
3/6
✓ Branch 1 taken 23484 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23484 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 23484 times.
✗ Branch 8 not taken.
 140904 arr.set_attribute("format", VTK::data_format_name(_encoder, _data_format));
70
2/4
✓ Branch 1 taken 23484 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23484 times.
✗ Branch 5 not taken.
 93936 arr.set_attribute("NumberOfComponents", ncomps);
71 46968 }
72
73 private:
74 const Encoder& _encoder;
75 const DataFormat& _data_format;
76 XMLElement& _element;
77 };
78
79 template<std::size_t dim>
80 class StructuredGridMapper {
81 public:
82 using Location = std::array<std::size_t, dim>;
83
84 4220 explicit StructuredGridMapper(std::vector<Location>&& map)
85 4220 : _map{std::move(map)}
86 4220 {}
87
88 35092 const auto& location(std::integral auto rank) const {
89 35092 return _map[rank];
90 }
91
92 35092 std::ranges::view auto ranks_below(const Location& loc, unsigned direction) const {
93 return std::views::filter(
94
1/2
✓ Branch 1 taken 17718 times.
✗ Branch 2 not taken.
 35092 std::views::iota(std::size_t{0}, _map.size()),
95 128124 [&, loc=loc, direction=direction] (const std::size_t rank) {
96 58660 const auto& rank_loc = _map[rank];
97
4/4
✓ Branch 0 taken 15862 times.
✓ Branch 1 taken 5110 times.
✓ Branch 2 taken 110744 times.
✓ Branch 3 taken 23140 times.
 154856 for (unsigned dir = 0; dir < dim; ++dir)
98
12/12
✓ Branch 0 taken 9064 times.
✓ Branch 1 taken 6798 times.
✓ Branch 4 taken 3954 times.
✓ Branch 5 taken 5110 times.
✓ Branch 6 taken 3954 times.
✓ Branch 7 taken 11908 times.
✓ Branch 8 taken 79340 times.
✓ Branch 9 taken 31404 times.
✓ Branch 12 taken 26456 times.
✓ Branch 13 taken 52884 times.
✓ Branch 14 taken 26456 times.
✓ Branch 15 taken 84288 times.
 126606 if (dir != direction && rank_loc[dir] != loc[dir]) {
99 30410 return false;
100 }
101 28250 return rank_loc[direction] < loc[direction];
102 }
103
1/2
✓ Branch 1 taken 17718 times.
✗ Branch 2 not taken.
 70184 );
104 }
105
106 private:
107 std::vector<Location> _map;
108 };
109
110 //! Helper for finding the locations of the sub-grids
111 //! associated with each rank in structured parallel grids
112 template<typename T, std::size_t dim>
113 class StructuredGridMapperHelper {
114 public:
115 using Origin = std::array<T, dim>;
116 using Basis = std::array<Origin, dim>;
117
118 4220 explicit StructuredGridMapperHelper(std::integral auto ranks,
119 const Basis& basis,
120 T default_epsilon = 1e-6)
121 4220 : _basis(basis)
122
1/2
✓ Branch 1 taken 2152 times.
✗ Branch 2 not taken.
 8440 , _origins(ranks)
123
1/2
✓ Branch 1 taken 2152 times.
✗ Branch 2 not taken.
 8440 , _set(ranks, false)
124 4220 , _default_epsilon{default_epsilon} {
125
1/2
✓ Branch 1 taken 2152 times.
✗ Branch 2 not taken.
 4220 std::ranges::fill(_reverse, false);
126 4220 }
127
128 5712 void reverse(int direction) {
129 5712 _reverse[direction] = !_reverse[direction];
130 5712 }
131
132 12588 void set_origin_for(int rank, Origin origin) {
133
2/4
✓ Branch 1 taken 6380 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 6380 times.
 12588 if (_set[rank])
134 ✗ throw ValueError("Origin for given rank already set");
135 12588 _origins[rank] = std::move(origin);
136
1/2
✓ Branch 1 taken 6380 times.
✗ Branch 2 not taken.
 12588 _set[rank] = true;
137 12588 }
138
139 4220 auto make_mapper() const {
140
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2152 times.
 4220 if (_origins.empty())
141 ✗ throw InvalidState("No origins have been set");
142
143
1/2
✓ Branch 2 taken 2152 times.
✗ Branch 3 not taken.
 4220 std::vector<std::array<std::size_t, dim>> map(_origins.size());
144
2/2
✓ Branch 1 taken 5956 times.
✓ Branch 2 taken 2152 times.
 15964 for (unsigned dir = 0; dir < dim; ++dir) {
145
1/2
✓ Branch 1 taken 5956 times.
✗ Branch 2 not taken.
 11744 auto ordinates = _get_ordinates(dir);
146
1/2
✓ Branch 1 taken 5956 times.
✗ Branch 2 not taken.
 11744 const auto eps = _epsilon(ordinates);
147
1/2
✓ Branch 2 taken 5956 times.
✗ Branch 3 not taken.
 11744 _sort_ordinates(ordinates, eps, _reverse[dir]);
148
149
2/2
✓ Branch 1 taken 3190 times.
✓ Branch 2 taken 2766 times.
 11744 if (ordinates.size() > 1) {
150
2/2
✓ Branch 1 taken 10532 times.
✓ Branch 2 taken 3190 times.
 27178 for (unsigned rank = 0; rank < _origins.size(); ++rank) {
151
1/2
✓ Branch 2 taken 10532 times.
✗ Branch 3 not taken.
 20884 const auto rank_origin_ordinate = _get_ordinate(_origins[rank], dir);
152
1/2
✓ Branch 1 taken 10532 times.
✗ Branch 2 not taken.
 36412 auto it = std::ranges::find_if(ordinates, [&] (const T o) {
153 using std::abs;
154 15798 return abs(o - rank_origin_ordinate) < eps;
155 });
156
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 10532 times.
 20884 if (it == ordinates.end())
157 ✗ throw InvalidState("Could not determine rank ordinate");
158 41768 map[rank][dir] = std::distance(ordinates.begin(), it);
159 }
160 } else {
161
2/2
✓ Branch 1 taken 7186 times.
✓ Branch 2 taken 2766 times.
 19658 for (unsigned rank = 0; rank < _origins.size(); ++rank)
162 14208 map[rank][dir] = 0;
163 }
164 }
165 8440 return StructuredGridMapper<dim>{std::move(map)};
166 4220 }
167
168 4220 Origin compute_origin() const {
169
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2152 times.
 4220 if (_origins.empty())
170 ✗ throw InvalidState("No origins have been set");
171
172 Origin result;
173
2/2
✓ Branch 0 taken 5956 times.
✓ Branch 1 taken 2152 times.
 15964 for (unsigned dir = 0; dir < dim; ++dir) {
174
1/2
✓ Branch 1 taken 5956 times.
✗ Branch 2 not taken.
 11744 auto ordinates = _get_ordinates(dir);
175
1/2
✓ Branch 2 taken 5956 times.
✗ Branch 3 not taken.
 11744 _sort_ordinates(ordinates, _reverse[dir]);
176 11744 result[dir] = ordinates[0];
177 }
178 4220 return result;
179 }
180
181 private:
182 23488 auto _get_ordinates(unsigned int axis) const {
183 23488 std::vector<T> result;
184
1/2
✓ Branch 2 taken 11912 times.
✗ Branch 3 not taken.
 23488 result.reserve(_origins.size());
185
2/4
✓ Branch 1 taken 11912 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11912 times.
✗ Branch 5 not taken.
 46976 std::ranges::copy(
186
1/2
✓ Branch 1 taken 11912 times.
✗ Branch 2 not taken.
 58236 std::views::transform(_origins, [&] (const Origin& o) {
187 35436 return _get_ordinate(o, axis);
188 }),
189 std::back_inserter(result)
190 );
191 23488 return result;
192 ✗ }
193
194 91068 auto _get_ordinate(const Origin& o, unsigned int axis) const {
195 91068 return dot_product(o, _basis[axis]);
196 }
197
198 11744 void _sort_ordinates(std::vector<T>& ordinates, bool reverse) const {
199
1/2
✓ Branch 1 taken 5956 times.
✗ Branch 2 not taken.
 11744 if (ordinates.size() > 1)
200
2/4
✓ Branch 1 taken 5956 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5956 times.
✗ Branch 5 not taken.
 11744 _sort_ordinates(ordinates, _epsilon(ordinates), reverse);
201 11744 }
202
203 23488 void _sort_ordinates(std::vector<T>& ordinates, const T& eps, bool reverse) const {
204
1/2
✓ Branch 1 taken 11912 times.
✗ Branch 2 not taken.
 23488 if (ordinates.size() > 1) {
205
3/6
✓ Branch 1 taken 11912 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11912 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 11912 times.
✗ Branch 8 not taken.
 46660 auto [it, _] = Ranges::sort_and_unique(ordinates, {}, [&] (T a, T b) {
206 using std::abs;
207 23524 return abs(a - b) < eps;
208 });
209
1/2
✓ Branch 4 taken 11912 times.
✗ Branch 5 not taken.
 23488 ordinates.erase(it, ordinates.end());
210
2/2
✓ Branch 0 taken 5712 times.
✓ Branch 1 taken 6200 times.
 23488 if (reverse)
211
1/2
✓ Branch 1 taken 5712 times.
✗ Branch 2 not taken.
 11424 std::ranges::reverse(ordinates);
212 }
213 23488 }
214
215 23488 T _epsilon(const std::vector<T>& ordinates) const {
216
1/2
✓ Branch 1 taken 11912 times.
✗ Branch 2 not taken.
 23488 const auto max = *std::ranges::max_element(ordinates);
217
1/2
✓ Branch 1 taken 11912 times.
✗ Branch 2 not taken.
 23488 const auto min = *std::ranges::min_element(ordinates);
218 23488 const auto size = max - min;
219
220
2/2
✓ Branch 0 taken 6380 times.
✓ Branch 1 taken 5532 times.
 23488 if (size > _default_epsilon) {
221
1/2
✓ Branch 2 taken 6380 times.
✗ Branch 3 not taken.
 12588 std::vector<T> dx; dx.reserve(ordinates.size());
222
2/4
✓ Branch 1 taken 6380 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 6380 times.
✗ Branch 7 not taken.
 12588 std::adjacent_difference(ordinates.begin(), ordinates.end(), std::back_inserter(dx));
223
1/2
✓ Branch 1 taken 6380 times.
✗ Branch 2 not taken.
 33652 std::ranges::for_each(dx, [] (auto& v) { v = std::abs(v); });
224
1/2
✓ Branch 1 taken 6380 times.
✗ Branch 2 not taken.
 12588 std::ranges::sort(dx);
225
1/2
✓ Branch 5 taken 16912 times.
✗ Branch 6 not taken.
 33472 for (T _dx : dx)
226
2/2
✓ Branch 0 taken 6380 times.
✓ Branch 1 taken 10532 times.
 33472 if (_dx > 1e-8*size)
227 12588 return _dx*0.1;
228
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 6380 times.
 12588 }
229
230 10900 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
241 template<typename Communicator>
242 class StructuredParallelGridHelper {
243 public:
244 6380 StructuredParallelGridHelper(const Communicator& comm, int root_rank = 0)
245 6380 : _comm{comm}
246 6380 , _root_rank{root_rank}
247 6380 {}
248
249 template<Concepts::Scalar CT,
250 std::size_t dim,
251 Concepts::StaticallySizedMDRange<2> B = std::array<std::array<CT, dim>, dim>>
252 12588 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
1/2
✓ Branch 1 taken 6380 times.
✗ Branch 2 not taken.
 12588 const auto num_ranks = Parallel::size(_comm);
257
1/2
✓ Branch 1 taken 6380 times.
✗ Branch 2 not taken.
 25176 std::vector<std::array<std::size_t, dim>> pieces_begin(num_ranks);
258
1/2
✓ Branch 1 taken 6380 times.
✗ Branch 2 not taken.
 12588 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
3/4
✓ Branch 1 taken 6380 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2152 times.
✓ Branch 4 taken 4228 times.
 12588 if (Parallel::rank(_comm) == _root_rank) {
263
3/6
✓ Branch 1 taken 2152 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2152 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2152 times.
✗ Branch 8 not taken.
 21594 const auto default_epsilon = 1e-6*std::ranges::max(all_origins | std::views::transform([] (CT v) {
264 using std::abs;
265 17718 return abs(v);
266 }));
267
1/2
✓ Branch 1 taken 2152 times.
✗ Branch 2 not taken.
 4220 const auto mapper_helper = _make_mapper_helper(basis, all_origins, is_negative_axis, default_epsilon);
268
1/2
✓ Branch 1 taken 2152 times.
✗ Branch 2 not taken.
 4220 const auto rank_mapper = mapper_helper.make_mapper();
269
1/2
✓ Branch 1 taken 2152 times.
✗ Branch 2 not taken.
 4220 origin = mapper_helper.compute_origin();
270
271
2/2
✓ Branch 0 taken 5956 times.
✓ Branch 1 taken 2152 times.
 15964 for (unsigned dir = 0; dir < dim; ++dir) {
272
2/2
✓ Branch 0 taken 17718 times.
✓ Branch 1 taken 5956 times.
 46836 for (int rank = 0; rank < num_ranks; ++rank) {
273
1/2
✓ Branch 2 taken 17718 times.
✗ Branch 3 not taken.
 35092 auto ranks_below = rank_mapper.ranks_below(rank_mapper.location(rank), dir);
274
3/6
✓ Branch 1 taken 17718 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 17718 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 17718 times.
✗ Branch 8 not taken.
 35092 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 5266 [&] (const std::size_t current, int r) {
279
2/4
✓ Branch 1 taken 1133 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4133 times.
✗ Branch 5 not taken.
 5266 return current + Parallel::access_gathered<dim>(all_extents, _comm, {dir, r});
280 }
281 );
282 35092 pieces_begin[rank][dir] = offset;
283
1/2
✓ Branch 1 taken 17718 times.
✗ Branch 2 not taken.
 35092 pieces_end[rank][dir] = offset + Parallel::access_gathered<dim>(all_extents, _comm, {dir, rank});
284 }
285
286
1/2
✓ Branch 3 taken 5956 times.
✗ Branch 4 not taken.
 23488 whole_extent[dir] = (*std::max_element(
287 pieces_end.begin(), pieces_end.end(),
288 11762 [&] (const auto& a1, const auto& a2) { return a1[dir] < a2[dir]; }
289 11744 ))[dir];
290 }
291 4220 }
292
293 return std::make_tuple(
294 12588 std::move(pieces_begin),
295 12588 std::move(pieces_end),
296 12588 std::move(whole_extent),
297 12588 std::move(origin)
298 25176 );
299 12588 }
300
301 private:
302 template<Concepts::StaticallySizedMDRange<2> Basis, Concepts::Scalar CT, std::size_t dim>
303 4220 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 4220 StructuredGridMapperHelper<CT, dim> helper(Parallel::size(_comm), basis, default_eps);
308
2/4
✓ Branch 1 taken 2152 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2152 times.
✗ Branch 5 not taken.
 16636 std::ranges::for_each(Parallel::ranks(_comm), [&] (int rank) {
309 6380 helper.set_origin_for(rank, Parallel::access_gathered<dim>(all_origins, _comm, rank));
310 });
311
312
2/2
✓ Branch 0 taken 5956 times.
✓ Branch 1 taken 2152 times.
 15964 for (unsigned dir = 0; dir < dim; ++dir)
313
2/2
✓ Branch 1 taken 2856 times.
✓ Branch 2 taken 3100 times.
 11744 if (is_negative_axis[dir])
314 5712 helper.reverse(dir);
315
316 4220 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_
326