dune-pdelab  2.5-dev
loadbalance.hh
Go to the documentation of this file.
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>
7 #include <dune/grid/common/datahandleif.hh>
8 #include <dune/grid/common/partitionset.hh>
9 
12 
13 namespace 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.
174  Buf bufWrapper(buff);
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.
229  Buf bufWrapper(buff);
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>
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>
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
static const int dim
Definition: adaptivity.hh:83
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
const Entity & e
Definition: localfunctionspace.hh:111
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
void loadBalanceGrid(Grid &grid, X &... x)
Load balance grid and restore gridfunctionspace and gridfunctions given as packs. ...
Definition: loadbalance.hh:499
Wrapper for message buffers of grid DataHandles that allows for sending different types of data...
Definition: polymorphicbufferwrapper.hh:26