GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/io/grid/vtkgriddatahandle.hh
Date: 2025-04-19 19:19:10
Exec Total Coverage
Lines: 108 223 48.4%
Functions: 6 63 9.5%
Branches: 127 483 26.3%
Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup InputOutput
10 * \brief A data handle for communicating grid data for VTK grids
11 */
12 #ifndef DUMUX_VTK_GRID_DATA_HANDLE_HH
13 #define DUMUX_VTK_GRID_DATA_HANDLE_HH
14
15 #include <memory>
16 #include <algorithm>
17 #include <map>
18 #include <numeric>
19 #include <array>
20 #include <vector>
21 #include <string>
22 #include <unordered_map>
23 #include <utility>
24 #include <iostream>
25 #include <ranges>
26
27 #include <dune/common/parallel/communication.hh>
28 #include <dune/common/parallel/future.hh>
29 #include <dune/geometry/dimension.hh>
30 #include <dune/grid/common/partitionset.hh>
31 #include <dune/grid/common/datahandleif.hh>
32
33 #include <dumux/io/vtk/vtkreader.hh>
34
35 namespace Dumux {
36
37 /*!
38 * \ingroup InputOutput
39 * \brief A data handle for communicating grid data for VTK grids
40 */
41 template<class Grid, class GridInput, class Data>
42 struct VtkGridDataHandle
43 : public Dune::CommDataHandleIF<VtkGridDataHandle<Grid, GridInput, Data>, typename Data::value_type>
44 {
45 using GridView = typename Grid::LevelGridView;
46
47
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 VtkGridDataHandle(const Grid& grid, const GridInput& gridInput, VTKReader::Data& cellData, VTKReader::Data& pointData)
48 4 : gridView_(grid.levelGridView(0))
49
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 , idSet_(grid.localIdSet())
50
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 , userCellData_(cellData)
51
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 , userPointData_(pointData)
52 {
53
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 const auto rank = grid.comm().rank();
54
55 // For the following to work we assume a sorted map of keys to values in the user data.
56 // This is not guaranteed by the VTKReader, so we need to sort the data first.
57
3/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 4 times.
 8 for (const auto& [key, data] : userCellData_)
58
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 cellData_[key] = std::move(userCellData_[key]);
59
3/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 4 times.
 8 for (const auto& [key, data] : userPointData_)
60
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 pointData_[key] = std::move(userPointData_[key]);
61
62 #if DUMUX_HAVE_GRIDFORMAT
63 // compute the number of components for each cell and point data array
64 // and the serialization size per entity
65 std::array<std::size_t, 2> numKeys{{ cellData_.size(), pointData_.size() }};
66 std::vector<std::size_t> keyComponents(numKeys[0] + numKeys[1], 0);
67 {
68 int n = 0;
69 for (const auto& [key, data] : cellData_)
70 keyComponents[n++] = rank == 0 ? data.size()/gridInput.numElements() : 0;
71 for (const auto& [key, data] : pointData_)
72 keyComponents[n++] = rank == 0 ? data.size()/gridInput.numVertices() : 0;
73 grid.comm().broadcast(keyComponents.data(), keyComponents.size(), 0);
74 }
75
76 const auto begin = keyComponents.begin();
77 cellDataComponents_.assign(begin, begin + numKeys[0]);
78 pointDataComponents_.assign(begin + numKeys[0], keyComponents.end());
79 numCellDataPerElement_ = std::accumulate(cellDataComponents_.begin(), cellDataComponents_.end(), 0UL);
80 numPointDataPerVertex_ = std::accumulate(pointDataComponents_.begin(), pointDataComponents_.end(), 0UL);
81 #else
82 // assume all data is on rank 0 (see grid manager)
83 // First broadcast how many keys we have
84
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 std::array<std::size_t, 2> numKeys{{ cellData_.size(), pointData_.size() }};
85
3/8
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✗ Branch 2 not taken.
✗ Branch 5 not taken.
 4 grid.comm().broadcast(numKeys.data(), 2, 0);
86
87 // then broadcast the length of the individual key strings
88 // and the number of component associated with each key (e.g. vector/tensor fields)
89
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 std::vector<std::size_t> keyLengthAndComponents(2*(numKeys[0] + numKeys[1]), 0);
90 {
91 4 int n = 0;
92
93 // key length
94
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 4 times.
 6 for (const auto& [key, data] : cellData_)
95 2 keyLengthAndComponents[n++] = key.size();
96
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 4 times.
 6 for (const auto& [key, data] : pointData_)
97 2 keyLengthAndComponents[n++] = key.size();
98
99 // number of components
100
3/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 4 times.
 6 for (const auto& [key, data] : cellData_)
101
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 keyLengthAndComponents[n++] = rank == 0 ? data.size()/gridInput.numElements() : 0;
102
3/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 4 times.
 6 for (const auto& [key, data] : pointData_)
103
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 keyLengthAndComponents[n++] = rank == 0 ? data.size()/gridInput.numVertices() : 0;
104
105
3/8
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
✗ Branch 2 not taken.
✗ Branch 5 not taken.
 4 grid.comm().broadcast(keyLengthAndComponents.data(), keyLengthAndComponents.size(), 0);
106 }
107
108 // save the number of components for each cell and point data array
109
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 const auto begin = keyLengthAndComponents.begin() + numKeys[0] + numKeys[1];
110
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 cellDataComponents_.assign(begin, begin + numKeys[0]);
111
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 pointDataComponents_.assign(begin + numKeys[0], keyLengthAndComponents.end());
112 4 numCellDataPerElement_ = std::accumulate(cellDataComponents_.begin(), cellDataComponents_.end(), 0UL);
113 8 numPointDataPerVertex_ = std::accumulate(pointDataComponents_.begin(), pointDataComponents_.end(), 0UL);
114
115 // then broadcast the actual keys
116
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 std::string keys; keys.resize(std::accumulate(keyLengthAndComponents.begin(), begin, 0UL));
117 {
118 4 int n = 0;
119
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 4 times.
 6 for (const auto& [key, data] : cellData_)
120
2/2
✓ Branch 0 taken 26 times.
✓ Branch 1 taken 2 times.
 28 for (const auto& c : key)
121 26 keys[n++] = c;
122
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 4 times.
 6 for (const auto& [key, data] : pointData_)
123
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 2 times.
 22 for (const auto& c : key)
124 20 keys[n++] = c;
125
126
2/5
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✗ Branch 2 not taken.
 4 grid.comm().broadcast(keys.data(), keys.size(), 0);
127 }
128
129 // create the entries in the cellData and pointData maps on all processes
130 ✗ std::size_t offset = 0;
131
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 8 for (int keyIdx = 0; keyIdx < numKeys[0]; ++keyIdx)
132 {
133
3/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
 4 if (std::string key{ keys, offset, keyLengthAndComponents[keyIdx] }; cellData_.count(key) == 0)
134
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 2 cellData_[key] = Data{};
135
136 4 offset += keyLengthAndComponents[keyIdx];
137 }
138
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 8 for (int keyIdx = numKeys[0]; keyIdx < numKeys[0] + numKeys[1]; ++keyIdx)
139 {
140
3/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
 4 if (std::string key{ keys, offset, keyLengthAndComponents[keyIdx] }; pointData_.count(key) == 0)
141
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 2 pointData_[key] = Data{};
142
143 4 offset += keyLengthAndComponents[keyIdx];
144 }
145 #endif
146
147 // write data into an id map
148
6/12
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 5088 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 5088 times.
✓ Branch 13 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 6 not taken.
 10192 for (const auto& element : elements(gridView_, Dune::Partitions::interior))
149 {
150
3/6
✓ Branch 1 taken 5088 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5088 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5088 times.
✗ Branch 8 not taken.
 5088 data_[idSet_.id(element)].resize(numCellDataPerElement_);
151
152 5088 int n = 0, l = 0;
153
2/2
✓ Branch 0 taken 5088 times.
✓ Branch 1 taken 5088 times.
 10176 for (const auto& [key, data] : cellData_)
154 {
155 5088 const auto nComp = cellDataComponents_[l++];
156
2/2
✓ Branch 0 taken 5088 times.
✓ Branch 1 taken 5088 times.
 10176 for (int k = 0; k < nComp; ++k)
157
4/8
✓ Branch 1 taken 5088 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5088 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5088 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 5088 times.
✗ Branch 11 not taken.
 5088 std::swap(cellData_[key][k + nComp*gridInput.insertionIndex(element)], data_[idSet_.id(element)][n++]);
158 }
159
160
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5088 times.
 5088 assert(n == numCellDataPerElement_);
161 }
162
163
6/12
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2650 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 4 times.
✓ Branch 13 taken 2650 times.
✗ Branch 3 not taken.
✗ Branch 6 not taken.
 5324 for (const auto& vertex : vertices(gridView_))
164 {
165
3/6
✓ Branch 1 taken 2650 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2650 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2650 times.
✗ Branch 8 not taken.
 2650 data_[idSet_.id(vertex)].resize(numPointDataPerVertex_);
166
167 2650 int n = 0, l = 0;
168
2/2
✓ Branch 0 taken 2650 times.
✓ Branch 1 taken 2650 times.
 5300 for (const auto& [key, data] : pointData_)
169 {
170 2650 const auto nComp = pointDataComponents_[l++];
171
2/2
✓ Branch 0 taken 2650 times.
✓ Branch 1 taken 2650 times.
 5300 for (int k = 0; k < nComp; ++k)
172
4/8
✓ Branch 1 taken 2650 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2650 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2650 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2650 times.
✗ Branch 11 not taken.
 2650 std::swap(pointData_[key][k + nComp*gridInput.insertionIndex(vertex)], data_[idSet_.id(vertex)][n++]);
173 }
174
175
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2650 times.
 2650 assert(n == numPointDataPerVertex_);
176 }
177 4 }
178
179 4 ~VtkGridDataHandle()
180 {
181 // resize arrays and unpack communicated data
182 4 const auto& indexSet = gridView_.indexSet();
183
184 {
185 4 int n = 0;
186
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 4 times.
 8 for (const auto& [key, data] : cellData_)
187 4 cellData_[key].resize(indexSet.size(0)*cellDataComponents_[n++]);
188 }
189 {
190 4 int n = 0;
191
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 4 times.
 8 for (const auto& [key, data] : pointData_)
192 4 pointData_[key].resize(indexSet.size(GridView::dimension)*pointDataComponents_[n++]);
193 }
194
195
2/5
✓ Branch 4 taken 5088 times.
✓ Branch 5 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 10188 for (const auto& element : elements(gridView_))
196 {
197 5088 int n = 0, l = 0;
198
2/2
✓ Branch 0 taken 5088 times.
✓ Branch 1 taken 5088 times.
 10176 for (const auto& [key, data] : cellData_)
199 {
200 5088 const auto nComp = cellDataComponents_[l++];
201
2/2
✓ Branch 0 taken 5088 times.
✓ Branch 1 taken 5088 times.
 10176 for (int k = 0; k < nComp; ++k)
202 5088 std::swap(cellData_[key][k + nComp*indexSet.index(element)], data_[idSet_.id(element)][n++]);
203 }
204 }
205
206
2/5
✓ Branch 4 taken 4 times.
✓ Branch 5 taken 2796 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 5604 for (const auto& vertex : vertices(gridView_))
207 {
208 2796 int n = 0, l = 0;
209
2/2
✓ Branch 0 taken 2796 times.
✓ Branch 1 taken 2796 times.
 5592 for (const auto& [key, data] : pointData_)
210 {
211 2796 const auto nComp = pointDataComponents_[l++];
212
2/2
✓ Branch 0 taken 2796 times.
✓ Branch 1 taken 2796 times.
 5592 for (int k = 0; k < nComp; ++k)
213 2796 std::swap(pointData_[key][k + nComp*indexSet.index(vertex)], data_[idSet_.id(vertex)][n++]);
214 }
215 }
216
217 // move data back from sorted internal storage to user data
218
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 4 times.
 8 for (const auto& [key, data] : cellData_)
219 4 userCellData_[key] = std::move(cellData_[key]);
220
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 4 times.
 8 for (const auto& [key, data] : pointData_)
221 4 userPointData_[key] = std::move(pointData_[key]);
222 4 }
223
224 Dune::CommDataHandleIF<VtkGridDataHandle<Grid, GridInput, Data>, typename Data::value_type>& interface()
225 { return *this; }
226
227 12 bool contains (int dim, int codim) const
228 12 { return codim == 0 || codim == dim; }
229
230 //! returns true if size per entity of given dim and codim is a constant
231 bool fixedSize(int dim, int codim) const
232 { return true; }
233
234 template<class Entity>
235 10468 std::size_t size (const Entity&) const
236 {
237 if constexpr (Entity::codimension == 0)
238
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 2544 return numCellDataPerElement_;
239 else if constexpr (Entity::codimension == GridView::dimension)
240
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 7924 return numPointDataPerVertex_;
241 else
242 return 0;
243 }
244
245 template<class MessageBufferImp, class Entity>
246 20936 void gather (MessageBufferImp& buff, const Entity& e) const
247 {
248 if constexpr (Entity::codimension == 0)
249 {
250 5088 const auto& data = data_[idSet_.id(e)];
251
2/2
✓ Branch 0 taken 2544 times.
✓ Branch 1 taken 2544 times.
 10176 for (int n = 0; n < numCellDataPerElement_; ++n)
252 5088 buff.write(data[n]);
253 }
254
255 if constexpr (Entity::codimension == GridView::dimension)
256 {
257 15848 const auto& data = data_[idSet_.id(e)];
258
2/2
✓ Branch 0 taken 7924 times.
✓ Branch 1 taken 7924 times.
 31696 for (int n = 0; n < numPointDataPerVertex_; ++n)
259 15848 buff.write(data[n]);
260 }
261 20936 }
262
263 template<class MessageBufferImp, class Entity>
264 20936 void scatter (MessageBufferImp& buff, const Entity& e, std::size_t n)
265 {
266 20936 auto& data = data_[idSet_.id(e)];
267 20936 data.resize(n);
268
269 if constexpr (Entity::codimension == 0)
270
2/2
✓ Branch 1 taken 2544 times.
✓ Branch 2 taken 2544 times.
 10176 for (int k = 0; k < numCellDataPerElement_; ++k)
271 5088 buff.read(data[k]);
272
273 if constexpr (Entity::codimension == GridView::dimension)
274
2/2
✓ Branch 1 taken 7924 times.
✓ Branch 2 taken 7924 times.
 31696 for (int k = 0; k < numPointDataPerVertex_; ++k)
275 15848 buff.read(data[k]);
276 20936 }
277
278 private:
279 using IdSet = typename Grid::LocalIdSet;
280
281 const GridView gridView_;
282 const IdSet &idSet_;
283
284 VTKReader::Data& userCellData_;
285 VTKReader::Data& userPointData_;
286
287 std::map<std::string, VTKReader::Data::mapped_type> cellData_;
288 std::map<std::string, VTKReader::Data::mapped_type> pointData_;
289
290 std::vector<std::size_t> cellDataComponents_;
291 std::vector<std::size_t> pointDataComponents_;
292
293 std::size_t numCellDataPerElement_;
294 std::size_t numPointDataPerVertex_;
295
296 mutable std::map< typename IdSet::IdType, std::vector<typename Data::value_type> > data_;
297 };
298
299 } // namespace Dumux
300
301 namespace Dumux::Detail::VtkData {
302
303 // for structured vtk data, we manually distribute the data to the other ranks
304 template<class Grid, class GridInput>
305 ✗ void communicateStructuredVtkData(const Grid& grid, const GridInput& gridInput, ::Dumux::VTKReader::Data& cellData, ::Dumux::VTKReader::Data& pointData)
306 {
307 #if HAVE_MPI // needed due to oversight in Dune::Communication interface
308 ✗ const auto commSize = grid.comm().size();
309 ✗ if (commSize <= 1)
310 ✗ return;
311
312 ✗ const auto rank = grid.comm().rank();
313
314 #ifndef NDEBUG
315 ✗ if (rank == 0)
316 {
317 ✗ std::cout << "Communicating structured VTK data...\n\n" << std::endl;
318 ✗ std::cout << "Grid has " << gridInput.numElements() << " elements and " << gridInput.numVertices() << " vertices." << std::endl;
319 }
320 #endif
321
322 // first some preliminary steps
323 // we need to sort the data
324 ✗ std::map<std::string, ::Dumux::VTKReader::Data::mapped_type> sortedCellData, sortedPointData;
325 ✗ for (const auto& [key, data] : cellData)
326 ✗ sortedCellData[key] = std::move(cellData[key]);
327 ✗ for (const auto& [key, data] : pointData)
328 ✗ sortedPointData[key] = std::move(pointData[key]);
329
330 // and we need to compute the number of components for each cell and point data array
331 // to know the message sizes
332 ✗ std::array<std::size_t, 2> numKeys{{ sortedCellData.size(), sortedPointData.size() }};
333 ✗ std::vector<std::size_t> keyComponents(numKeys[0] + numKeys[1], 0);
334 {
335 ✗ int n = 0;
336 ✗ for (const auto& [key, data] : sortedCellData)
337 ✗ keyComponents[n++] = rank == 0 ? data.size()/gridInput.numElements() : 0;
338 ✗ for (const auto& [key, data] : sortedPointData)
339 ✗ keyComponents[n++] = rank == 0 ? data.size()/gridInput.numVertices() : 0;
340 ✗ grid.comm().broadcast(keyComponents.data(), keyComponents.size(), 0);
341 }
342
343 ✗ const auto begin = keyComponents.begin();
344 ✗ const std::size_t numCellDataPerElement = std::accumulate(begin, begin + numKeys[0], 0UL);
345 ✗ const std::size_t numPointDataPerVertex = std::accumulate(begin + numKeys[0], keyComponents.end(), 0UL);
346
347 #ifndef NDEBUG
348 ✗ std::cout << "Rank " << rank << ": numbers of components for each cell and point data array: "
349 ✗ << numCellDataPerElement << " " << numPointDataPerVertex << std::endl;
350 #endif
351
352 // each process knows:
353 // - total number of cells and vertices
354 // - its data indices for each cell and vertex (global numbering)
355 // process 0 has all the data
356
357 // each rank decides which data they need
358 ✗ const auto& gridView = grid.levelGridView(0);
359 ✗ std::vector<std::size_t> requestedData(gridView.size(0) + gridView.size(Grid::dimension));
360 ✗ const std::size_t numElements = gridView.size(0);
361 ✗ const std::size_t numVertices = gridView.size(Grid::dimension);
362 ✗ for (const auto& element : elements(gridView))
363 {
364 ✗ requestedData[gridView.indexSet().index(element)] = gridInput.insertionIndex(element);
365 ✗ for (const auto& vertex : subEntities(element, Dune::Codim<Grid::dimension>{}))
366 ✗ requestedData[numElements + gridView.indexSet().index(vertex)] = gridInput.insertionIndex(vertex);
367 }
368
369 // gather the sizes of the data on rank 0
370 ✗ std::vector<std::size_t> numData_;
371 ✗ if (rank == 0)
372 ✗ numData_.resize(2*commSize);
373 ✗ std::array<std::size_t, 2> localNumData{{ numElements, numVertices }};
374 ✗ grid.comm().gather(localNumData.data(), numData_.data(), 2, 0);
375
376 #ifndef NDEBUG
377 ✗ if (rank == 0)
378 {
379 ✗ std::cout << "Number of elements and vertices on each rank: ";
380 ✗ for (std::size_t i = 0; i < commSize; ++i)
381 ✗ std::cout << numData_[2*i] << " " << numData_[2*i + 1] << " ";
382 ✗ std::cout << std::endl;
383 }
384 #endif
385
386 // send the data request to rank 0
387 using RequestData = std::vector<std::size_t>;
388 using FutureIndices = Dune::Future<RequestData>;
389 ✗ std::unique_ptr<FutureIndices> sendRequest;
390 if (rank != 0)
391 ✗ sendRequest = std::make_unique<FutureIndices>(
392 ✗ std::move(grid.comm().isend(std::move(requestedData), 0, /*tag*/0))
393 );
394
395 // receive the data on rank 0
396 ✗ std::vector<RequestData> requestedDataAll(commSize);
397 ✗ if (rank == 0)
398 {
399 ✗ std::vector<std::unique_ptr<FutureIndices>> receiveRequests(commSize-1);
400 ✗ for (std::size_t i = 0; i < commSize; ++i)
401 {
402 ✗ requestedDataAll[i].resize(numData_[2*i] + numData_[2*i + 1]);
403
404 ✗ if (i == 0)
405 ✗ requestedDataAll[i] = std::move(requestedData);
406 else
407 {
408 ✗ receiveRequests[i-1] = std::make_unique<FutureIndices>(
409 ✗ std::move(grid.comm().irecv(requestedDataAll[i], i, /*tag*/0))
410 );
411 }
412 }
413
414 /// TODO: actually we want to call MPI_Waitall here: how to do this with the futures?
415 ✗ std::ranges::for_each(receiveRequests, [](auto& request) { request->wait(); });
416
417 #ifndef NDEBUG
418 ✗ std::cout << "Received data from all ranks." << std::endl;
419 ✗ for (std::size_t i = 0; i < commSize; ++i)
420 ✗ std::cout << "From rank " << i << ": " << requestedDataAll[i].size() << " index size." << std::endl;
421 #endif
422 ✗ }
423
424 // Now we need to pack the data and send it to the other ranks
425 using PackedData = std::vector<double>;
426 ✗ PackedData receivedFlatData(localNumData[0]*numCellDataPerElement + localNumData[1]*numPointDataPerVertex);
427 ✗ if (rank == 0)
428 {
429 using FutureData = Dune::Future<PackedData>;
430 ✗ std::vector<std::unique_ptr<FutureData>> sendRequests(commSize-1);
431 ✗ for (std::size_t i = 0; i < commSize; ++i)
432 {
433 ✗ const auto& requestedIndices = requestedDataAll[i];
434 ✗ PackedData packedData(numData_[2*i]*numCellDataPerElement + numData_[2*i + 1]*numPointDataPerVertex);
435
436 // pack the data
437 ✗ std::size_t n = 0, l = 0;
438 ✗ for (const auto& [key, data] : sortedCellData)
439 {
440 ✗ const auto nComp = keyComponents[l++];
441 ✗ for (std::size_t k = 0; k < numData_[2*i]; ++k)
442 {
443 ✗ const auto idx = requestedIndices[k];
444 ✗ for (std::size_t j = 0; j < nComp; ++j)
445 ✗ packedData[n++] = data[j + nComp*idx];
446 }
447 }
448
449 ✗ const auto pointDataOffsett = numData_[2*i];
450 ✗ for (const auto& [key, data] : sortedPointData)
451 {
452 ✗ const auto nComp = keyComponents[l++];
453 ✗ for (std::size_t k = 0; k < numData_[2*i + 1]; ++k)
454 {
455 ✗ const auto idx = requestedIndices[k + pointDataOffsett];
456 ✗ for (std::size_t j = 0; j < nComp; ++j)
457 ✗ packedData[n++] = data[j + nComp*idx];
458 }
459 }
460
461 ✗ if (n != packedData.size())
462 ✗ DUNE_THROW(Dune::InvalidStateException, "Packed data size does not match expected size.");
463
464 ✗ if (i == 0)
465 ✗ receivedFlatData = std::move(packedData);
466 else
467 {
468 // send the data to rank i
469 ✗ sendRequests[i-1] = std::make_unique<FutureData>(
470 ✗ std::move(grid.comm().isend(std::move(packedData), i, /*tag*/1))
471 );
472 }
473 }
474
475 /// TODO: actually we want to call MPI_Waitall here: how to do this with the futures?
476 ✗ std::ranges::for_each(sendRequests, [](auto& request) { request->wait(); });
477 ✗ }
478 else
479 {
480 // receive the data from rank 0
481 ✗ auto receiveRequest = grid.comm().irecv(receivedFlatData, 0, /*tag*/1);
482
483 // finalize all communication
484 ✗ sendRequest->wait();
485 ✗ receiveRequest.wait();
486 ✗ }
487
488 #ifndef NDEBUG
489 ✗ std::cout << "On rank " << rank << ", the received data size is " << receivedFlatData.size() << std::endl;
490 #endif
491
492 // unpack the data on each process into cellData and pointData
493 ✗ std::size_t n = 0, l = 0;
494 ✗ for (const auto& [key, data] : sortedCellData)
495 {
496 ✗ auto& out = cellData[key];
497 ✗ const auto nComp = keyComponents[l++];
498 ✗ out.resize(localNumData[0]*nComp);
499 ✗ for (std::size_t k = 0; k < localNumData[0]; ++k)
500 ✗ for (std::size_t j = 0; j < nComp; ++j)
501 ✗ out[j + nComp*k] = receivedFlatData[n++];
502 }
503
504 ✗ for (const auto& [key, data] : sortedPointData)
505 {
506 ✗ auto& out = pointData[key];
507 ✗ const auto nComp = keyComponents[l++];
508 ✗ out.resize(localNumData[1]*nComp);
509 ✗ for (std::size_t k = 0; k < localNumData[1]; ++k)
510 ✗ for (std::size_t j = 0; j < nComp; ++j)
511 ✗ out[j + nComp*k] = receivedFlatData[n++];
512 }
513
514 ✗ if (n != receivedFlatData.size())
515 ✗ DUNE_THROW(Dune::InvalidStateException, "Unpacked data size does not match expected size.");
516
517 #ifndef NDEBUG
518 ✗ if (rank == 0)
519 ✗ std::cout << "\n\n+++++++++++++++++++++++++++++++++++++++++++++" << std::endl;
520 #endif
521 #endif
522 ✗ }
523
524 } // end namespace Dumux::Detail::VtkData
525
526 #endif
527