DUNE PDELab (git)

loadbalance.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_PDELAB_GRIDFUNCTIONSPACE_LOADBALANCE_HH
4#define DUNE_PDELAB_GRIDFUNCTIONSPACE_LOADBALANCE_HH
5
6#include <dune/geometry/dimension.hh>
8#include <dune/grid/common/partitionset.hh>
9
10#include<dune/pdelab/common/polymorphicbufferwrapper.hh>
11#include <dune/pdelab/gridfunctionspace/entityindexcache.hh>
12
13namespace Dune {
14 namespace PDELab {
15
16#ifndef DOXYGEN
17 namespace impl{
18
23 template<typename... T>
24 class LoadBalanceDataHandle
25 : public Dune::CommDataHandleIF<LoadBalanceDataHandle<T...>,
26 typename std::decay<typename std::tuple_element<0,std::tuple<T...>>::type>::type
27 ::Map::mapped_type::value_type>
28 {
29 public:
30
31 using R0 = typename std::decay<typename std::tuple_element<0,std::tuple<T...>>::type>::type
32 ::Map::mapped_type::value_type;
33
34 LoadBalanceDataHandle(std::tuple<T&...>&& mapTuple)
35 : _mapTuple(mapTuple)
36 {
37 }
38
40 bool contains (int dim, int codim) const
41 {
42 return true;
43 }
44
46 bool fixedSize (int dim, int codim) const
47 {
48 // We return false here, since gather and scatter is called
49 // for all entities of all levels of the grid but we only
50 // communicate data for a given gridView. This means there
51 // will be gather and scatter calls for cells that don't exist
52 // in our current grid view and we won't communicate anything
53 // for these cells.
54 return false;
55 }
56
57 // End of tmp recursion
58 template <std::size_t I, typename EntityType>
59 inline typename std::enable_if<I==sizeof...(T),void>::type
60 sizeTMP(std::size_t& commSize, EntityType& e) const
61 {
62 }
63
64 // Go through tuple and accumulate communication size
65 template <std::size_t I=0, typename EntityType>
66 inline typename std::enable_if<I<sizeof...(T),void>::type
67 sizeTMP(std::size_t& commSize, EntityType& e) const
68 {
69 // Compare if data type for this entry equals data typo of first entry
70 using R = typename std::decay<typename std::tuple_element<I,std::tuple<T...>>::type>::type
71 ::Map::mapped_type::value_type;
72 static_assert(std::is_same<R,R0>::value,"Different field type of vectors not supported by DH");
73
74 // Current function space
75 auto& gfs = std::get<I>(_mapTuple)._gfs;
76
77 // Only communicate if there are dofs for this codimension.
78 if (gfs.finiteElementMap().hasDOFs(e.codimension)){
79
80 // We first have to communicate the size of data we send for this particular vector
81 commSize += sizeof(R);
82
83 // Only communicate data if entity lies in our entitySet
84 if (gfs.entitySet().contains(e)){
85 // Get some types
86 using GFS = typename std::decay<decltype(gfs)>::type;
87 using EntitySet = typename GFS::Traits::EntitySet;
88 using IDSet = typename EntitySet::Traits::GridView::Grid::LocalIdSet;
89
90 // Find element id in map
91 const IDSet& idSet(gfs.entitySet().gridView().grid().localIdSet());
92 auto& map = std::get<I>(_mapTuple)._map;
93 auto find = map.find(idSet.id(e));
94 assert (find!=map.end());
95
96 // Add size of degrees of freedom vetor
97 commSize += find->second.size()*sizeof(R);
98 }
99 }
100
101 // Get size for next vector
102 sizeTMP<I+1>(commSize,e);
103 }
104
106 template<class EntityType>
107 std::size_t size (EntityType& e) const
108 {
109 // Use tmp to sum up sizes for all vectors
110 std::size_t commSize(0.0);
111 sizeTMP<0>(commSize,e);
112 return commSize;
113 }
114
115
116 // End of tmp recursion
117 template <std::size_t I, typename Buf, typename Entity>
118 inline typename std::enable_if<I==sizeof...(T),void>::type
119 gatherTMP(Buf& buf, const Entity& e) const
120 {
121 }
122
123 // Go through tuple and call gather for all vectors
124 template <std::size_t I=0, typename Buf, typename Entity>
125 inline typename std::enable_if<I<sizeof...(T),void>::type
126 gatherTMP(Buf& buf, const Entity& e) const
127 {
128 // Current gridFunctionSpace and dof map for the vector
129 auto& gfs = std::get<I>(_mapTuple)._gfs;
130 auto& map = std::get<I>(_mapTuple)._map;
131
132 // Communication is called for every level and every entity of
133 // every codim of all entities where load balance will result in
134 // changes. We only send dofs for current gridView/entitySet.
135 if (gfs.finiteElementMap().hasDOFs(e.codimension)){
136 if (gfs.entitySet().contains(e)){
137 // Get important types
138 using GFS = typename std::decay<decltype(gfs)>::type;
139 using EntitySet = typename GFS::Traits::EntitySet;
140 using IDSet = typename EntitySet::Traits::GridView::Grid::LocalIdSet;
141
142 // Find element id in map
143 const IDSet& idSet(gfs.entitySet().gridView().grid().localIdSet());
144
145 // Find entity id in map.
146 auto find = map.find(idSet.id(e));
147 assert (find!=map.end());
148
149 // Send size we need to communicate for this vector
150 buf.write (static_cast<R0>(find->second.size()));
151
152 // Send dofs
153 for (size_t i=0; i<find->second.size(); ++i){
154 buf.write(find->second[i]);
155 }
156 }
157 else {
158 // Only communicate that we don't communicate any DOFs
159 R0 tmp(0);
160 buf.write (tmp);
161 }
162 } // hasDOFs
163
164 // Call gather for next vector
165 gatherTMP<I+1> (buf,e);
166 }
167
169 template<class MessageBuffer, class EntityType>
170 void gather (MessageBuffer& buff, const EntityType& e) const
171 {
172 // Communicate different things than char.
175
176 // Call gather for all vectors using tmp
177 gatherTMP<0> (bufWrapper,e);
178 }
179
180 // End of tmp reucrsion
181 template <std::size_t I, typename Buf, typename Entity>
182 inline typename std::enable_if<I==sizeof...(T),void>::type
183 scatterTMP(Buf& buf, const Entity& e) const
184 {
185 }
186
187 // Go through tuple receive DOFs
188 template <std::size_t I=0, typename Buf, typename Entity>
189 inline typename std::enable_if<I<sizeof...(T),void>::type
190 scatterTMP(Buf& buf, const Entity& e) const
191 {
192 auto& gfs = std::get<I>(_mapTuple)._gfs;
193 auto& map = std::get<I>(_mapTuple)._map;
194 if (gfs.finiteElementMap().hasDOFs(e.codimension)){
195
196 // Receive number of DOFs for this vector
197 R0 tmp;
198 buf.read(tmp);
199 std::size_t numberOfEntries(0);
200 numberOfEntries = (size_t) tmp;
201
202 // Create vector of DOFs and receive DOFs
203 std::vector<R0> dofs(numberOfEntries);
204 for (size_t i=0; i<numberOfEntries; ++i){
205 buf.read(dofs[i]);
206 }
207
208 // Store id and dofs in map
209 const auto& id_set = gfs.entitySet().grid().localIdSet();
210 map.insert({{id_set.id(e),dofs}});
211 }
212
213 // Call scatter for next vetcor using tmp
214 scatterTMP<I+1> (buf,e);
215 }
216
217
224 template<class MessageBuffer, class EntityType>
225 void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
226 {
227 // Communicate different things than char.
230
231 // Call scatter for all vectors using tmp
232 scatterTMP<0> (bufWrapper, e);
233 }
234
235 private:
236 // Tuple storing GFSandMap structs for every vector that should get adapted
237 std::tuple<T&...>& _mapTuple;
238 };
239
240
242 template <typename GFS, typename V, typename MAP, int codim>
243 void loadBalanceMapFiller (const GFS& gfs, V& v, MAP& map)
244 {
245 using IndexCache = Dune::PDELab::EntityIndexCache<GFS>;
246 using LocalView = typename V::template LocalView<IndexCache>;
247 IndexCache indexCache(gfs);
248 LocalView localView(v);
249 const auto& id_set = gfs.entitySet().grid().localIdSet();
250
251 // Iterate over all interiorBorder entities of codim.
252 for (const auto& e : entities(gfs.entitySet(),Dune::Codim<codim>(),Dune::Partitions::interiorBorder)){
253 // Bind cache to entity.
254 indexCache.update(e);
255 localView.bind(indexCache);
256
257 // Store dofs of entity in std::vector.
258 std::vector<typename LocalView::ElementType> dofs;
259 for (std::size_t i=0; i<localView.size(); ++i){
260 dofs.push_back(localView[i]);
261 }
262
263 // Insert id and vector in map.
264 map.insert ( {{id_set.id(e),dofs}});
265
266 // Unbind cache.
267 localView.unbind();
268 }
269 }
270
272 template <int codim>
273 struct FillLoadBalanceDOFMap
274 {
275 template <typename GFS, typename V, typename MAP>
276 static void fillMap (const GFS& gfs, V& v, MAP& map)
277 {
278 if (gfs.finiteElementMap().hasDOFs(codim)){
279 loadBalanceMapFiller<GFS,V,MAP,codim>(gfs,v,map);
280 }
281 FillLoadBalanceDOFMap<codim-1>::fillMap(gfs,v,map);
282 }
283 };
284 template <>
285 struct FillLoadBalanceDOFMap<0>
286 {
287 template <typename GFS, typename V, typename MAP>
288 static void fillMap (const GFS& gfs, V& v, MAP& map)
289 {
290 if (gfs.finiteElementMap().hasDOFs(0)){
291 loadBalanceMapFiller<GFS,V,MAP,0>(gfs,v,map);
292 }
293 }
294 };
295
297 template <typename GFS, typename V, typename MAP, int codim>
298 void loadBalanceMapReader (const GFS& gfs, V& v, MAP& map)
299 {
300 using IndexCache = Dune::PDELab::EntityIndexCache<GFS>;
301 using LocalView = typename V::template LocalView<IndexCache>;
302 IndexCache indexCache(gfs);
303 LocalView localView(v);
304 const auto& id_set = gfs.entitySet().grid().localIdSet();
305
306 // Iterate over all interiorBorder entities of codim.
307 for (const auto& e : entities(gfs.entitySet(),Dune::Codim<codim>(),Dune::Partitions::interiorBorder)){
308 // Bind cache to entity.
309 indexCache.update(e);
310 localView.bind(indexCache);
311
312 // Find key in map and get vector of dofs.
313 auto find = map.find(id_set.id(e));
314 auto& dofs(find->second);
315
316 // Assert that we found element and that sizes of dof vector and local view match.
317 assert(find!=map.end());
318 assert(dofs.size()==localView.size());
319
320 // Store dofs in local view.
321 for (std::size_t i=0; i<dofs.size(); ++i){
322 localView[i]=dofs[i];
323 }
324
325 // Write changes to underlying container.
326 localView.commit();
327
328 // Unbind cache.
329 localView.unbind();
330 }
331 }
332
334 template <int codim>
335 struct ReadLoadBalanceDOFMap
336 {
337 template <typename GFS, typename V, typename MAP>
338 static void readMap (const GFS& gfs, V& v, MAP& map)
339 {
340 if (gfs.finiteElementMap().hasDOFs(codim)){
341 loadBalanceMapReader<GFS,V,MAP,codim>(gfs,v,map);
342 }
343 ReadLoadBalanceDOFMap<codim-1>::readMap(gfs,v,map);
344 }
345 };
346 template <>
347 struct ReadLoadBalanceDOFMap<0>
348 {
349 template <typename GFS, typename V, typename MAP>
350 static void readMap (const GFS& gfs, V& v, MAP& map)
351 {
352 if (gfs.finiteElementMap().hasDOFs(0)){
353 loadBalanceMapReader<GFS,V,MAP,0>(gfs,v,map);
354 }
355 }
356 };
357
358
359 // Store reference to function space and map
360 template <typename G, typename M>
361 struct GFSAndMap
362 {
363 // Export types
364 using GFS = G;
365 using Map = M;
366
367 GFSAndMap (GFS& gfs, Map& m) : _gfs(gfs), _map(m)
368 {
369 }
370
371 GFS& _gfs;
372 Map& _map;
373 };
374
375 // Create a GFSAndMap struct
376 template <typename GFS, typename M>
377 GFSAndMap<GFS,M> packGFSAndMap(GFS& gfs, M& m)
378 {
379 GFSAndMap<GFS,M> pack(gfs,m);
380 return pack;
381 }
382
383
384 // Forward declarations needed for the tmp recursion
385 template <typename Grid, typename... T>
386 void iteratePacks(Grid& grid, std::tuple<T&...>&& mapTuple);
387 template <typename Grid, typename... T, typename X, typename... XS>
388 void iteratePacks(Grid& grid, std::tuple<T&...>&& mapTuple, X& x, XS&... xs);
389
390 // This function is called after the last vector of the tuple. Here
391 // the next pack is called. On the way back we update the current
392 // function space.
393 template<std::size_t I = 0, typename Grid, typename... T, typename X, typename... XS>
394 inline typename std::enable_if<I == std::tuple_size<typename X::Tuple>::value, void>::type
395 iterateTuple(Grid& grid, std::tuple<T&...>&& mapTuple, X& x, XS&... xs)
396 {
397 // Iterate next pack
398 iteratePacks(grid,std::move(mapTuple),xs...);
399
400 // On our way back we need to update the current function space
401 x._gfs.update(true);
402 }
403
404 /* In this function we store the data of the current vector (indicated
405 * by template parameter I) of the current pack. After recursively
406 * iterating through the other packs and vectors we replay the data.
407 *
408 * @tparam I std:size_t used for tmp
409 * @tparam Grid Grid type
410 * @tparam T... Types of tuple elements (for storing data transfer maps)
411 * @tparam X Current pack
412 * @tparam ...XS Remaining packs
413 */
414 template<std::size_t I = 0, typename Grid, typename... T, typename X, typename... XS>
415 inline typename std::enable_if<I < std::tuple_size<typename X::Tuple>::value, void>::type
416 iterateTuple(Grid& grid, std::tuple<T&...>&& mapTuple, X& x, XS&... xs)
417 {
418 // Get some basic types
419 using GFS = typename X::GFS;
420 using Tuple = typename X::Tuple;
421 using V = typename std::decay<typename std::tuple_element<I,Tuple>::type>::type;
422 using IDSet = typename Grid::LocalIdSet;
423 using ID = typename IDSet::IdType;
424 using R = typename V::field_type;
425
426 // Store vector of dofs for every entitiy of every codim.
427 using MAP = std::unordered_map<ID,std::vector<R>>;
428 MAP map;
429 FillLoadBalanceDOFMap<GFS::Traits::GridView::dimension>::fillMap(x._gfs,std::get<I>(x._tuple),map);
430
431 // Pack gfs and tuple
432 auto mapPack = packGFSAndMap(x._gfs,map);
433 auto newMapTuple = std::tuple_cat(mapTuple,std::tie(mapPack));
434
435 // Recursively iterate through remaining vectors (and packs). Load
436 // balancing will be done at the end of recursion.
437 iterateTuple<I+1>(grid,std::move(newMapTuple),x,xs...);
438
439 // Restore solution from dof map.
440 std::get<I>(x._tuple) = V(x._gfs,0.0);
441 ReadLoadBalanceDOFMap<GFS::Traits::GridView::dimension>::readMap(x._gfs,std::get<I>(x._tuple),map);
442 }
443
444 template <typename... T>
445 LoadBalanceDataHandle<T...> createLoadBalanceDataHandle (std::tuple<T&...>&& mapTuple)
446 {
447 LoadBalanceDataHandle<T...> dh(std::move(mapTuple));
448 return dh;
449 }
450
451 // This gets called after the last pack. After this function call we
452 // have visited every vector of every pack and we will go back through
453 // the recursive function calls.
454 template <typename Grid, typename... T>
455 void iteratePacks(Grid& grid, std::tuple<T&...>&& mapTuple)
456 {
457 // Create data handle and use load balancing with communication.
458 auto dh = createLoadBalanceDataHandle(std::move(mapTuple));
459
460 std::cout << "Calling load balance with data communication" << std::endl;
461 grid.loadBalance(dh);
462 }
463
464 /* Use template meta programming to iterate over packs at compile time
465 *
466 * In order to adapt our grid and all vectors of all packs we need to
467 * do the following:
468 * - Iterate over all vectors of all packs.
469 * - Store the data from the vectors where things could change.
470 * - Load Balance our grid.
471 * - Update function spaces and restore data.
472 *
473 * The key point is that we need the object that stores the data to
474 * replay it. Because of that we can not just iterate over the packs
475 * and within each pack iterate over the vectors but we have to make
476 * one big recursion. Therefore we iterate over the vectors of the
477 * current pack.
478 */
479 template <typename Grid, typename... T, typename X, typename... XS>
480 void iteratePacks(Grid& grid, std::tuple<T&...>&& mapTuple, X& x, XS&... xs)
481 {
482 iterateTuple(grid,std::move(mapTuple),x,xs...);
483 }
484
485 } // namespace impl
486#endif // DOXYGEN
487
488
498 template <typename Grid, typename... X>
499 void loadBalanceGrid(Grid& grid, X&... x)
500 {
501 // Create tuple where all data transfer maps get stored
502 std::tuple<> mapTuple;
503
504 // Iterate over packs
505 impl::iteratePacks(grid,std::move(mapTuple),x...);
506 }
507
508 } // namespace PDELab
509} // namespace Dune
510
511#endif // DUNE_PDELAB_GRIDFUNCTIONSPACE_LOADBALANCE_HH
CommDataHandleIF describes the features of a data handle for communication in parallel runs using the...
Definition: datahandleif.hh:78
GridFamily::Traits::LocalIdSet LocalIdSet
A type that is a model of Dune::IdSet which provides a unique and persistent numbering for all entiti...
Definition: grid.hh:509
Wrapper for message buffers of grid DataHandles that allows for sending different types of data.
Definition: polymorphicbufferwrapper.hh:32
Describes the parallel communication interface class for MessageBuffers and DataHandles.
constexpr InteriorBorder interiorBorder
PartitionSet for the interior and border partitions.
Definition: partitionset.hh:286
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
constexpr std::bool_constant<((II==value)||...)> contains(std::integer_sequence< T, II... >, std::integral_constant< T, value >)
Checks whether or not a given sequence contains a value.
Definition: integersequence.hh:137
Static tag representing a codimension.
Definition: dimension.hh:24
std::size_t fixedSize
The number of data items per index if it is fixed, 0 otherwise.
Definition: variablesizecommunicator.hh:264
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)