1/*
2 * SpaceDistribution.hpp
3 *
4 * Created on: Dec 3, 2016
5 * Author: i-bird
6 */
7
8#ifndef SRC_DECOMPOSITION_DISTRIBUTION_SPACEDISTRIBUTION_HPP_
9#define SRC_DECOMPOSITION_DISTRIBUTION_SPACEDISTRIBUTION_HPP_
10
11#include "util/mathutil.hpp"
12#include "NN/CellList/CellDecomposer.hpp"
13#include "Grid/grid_key_dx_iterator_hilbert.hpp"
14
15/*! \brief Class that distribute sub-sub-domains across processors using an hilbert curve
16 * to divide the space
17 *
18 * ### Initialize a Cartesian graph and decompose
19 * \snippet Distribution_unit_tests.hpp Initialize a Space Cartesian graph and decompose
20 *
21 *
22 */
23template<unsigned int dim, typename T>
24class SpaceDistribution
25{
26 //! Vcluster
27 Vcluster<> & v_cl;
28
29 //! Structure that store the cartesian grid information
30 grid_sm<dim, void> gr;
31
32 //! rectangular domain to decompose
33 Box<dim, T> domain;
34
35 //! Global sub-sub-domain graph
36 Graph_CSR<nm_v<dim>, nm_e> gp;
37
38
39public:
40
41 /*! Constructor
42 *
43 * \param v_cl Vcluster to use as communication object in this class
44 */
45 SpaceDistribution(Vcluster<> & v_cl)
46 :v_cl(v_cl)
47 {
48 }
49
50 /*! Copy constructor
51 *
52 * \param pm Distribution to copy
53 *
54 */
55 SpaceDistribution(const ParMetisDistribution<dim,T> & pm)
56 :v_cl(pm.v_cl)
57 {
58 this->operator=(pm);
59 }
60
61 /*! Copy constructor
62 *
63 * \param pm Distribution to copy
64 *
65 */
66 SpaceDistribution(SpaceDistribution<dim,T> && pm)
67 :v_cl(pm.v_cl)
68 {
69 this->operator=(pm);
70 }
71
72 /*! \brief Create the Cartesian graph
73 *
74 * \param grid info
75 * \param dom domain
76 */
77 void createCartGraph(grid_sm<dim, void> & grid, Box<dim, T> dom)
78 {
79 size_t bc[dim];
80
81 for (size_t i = 0 ; i < dim ; i++)
82 bc[i] = NON_PERIODIC;
83
84 // Set grid and domain
85 gr = grid;
86 domain = dom;
87
88 // Create a cartesian grid graph
89 CartesianGraphFactory<dim, Graph_CSR<nm_v<dim>, nm_e>> g_factory_part;
90 gp = g_factory_part.template construct<NO_EDGE, nm_v_id, T, dim - 1, 0>(gr.getSize(), domain, bc);
91
92 // Init to 0.0 axis z (to fix in graphFactory)
93 if (dim < 3)
94 {
95 for (size_t i = 0; i < gp.getNVertex(); i++)
96 gp.vertex(i).template get<nm_v_x>()[2] = 0.0;
97 }
98 for (size_t i = 0; i < gp.getNVertex(); i++)
99 gp.vertex(i).template get<nm_v_global_id>() = i;
100
101 }
102
103 /*! \brief Get the current graph (main)
104 *
105 */
106 Graph_CSR<nm_v<dim>, nm_e> & getGraph()
107 {
108 return gp;
109 }
110
111 /*! \brief Create the decomposition
112 *
113 */
114 void decompose()
115 {
116 // Get the number of processing units
117 size_t Np = v_cl.getProcessingUnits();
118
119 // Calculate the best number of sub-domains for each
120 // processor
121 size_t N_tot = gr.size();
122 size_t N_best_each = N_tot / Np;
123 size_t N_rest = N_tot % Np;
124
125 openfpm::vector<size_t> accu(Np);
126 accu.get(0) = N_best_each + ((0 < N_rest)?1:0);
127 for (size_t i = 1 ; i < Np ; i++)
128 accu.get(i) = accu.get(i-1) + N_best_each + ((i < N_rest)?1:0);
129
130 // Get the maximum along dimensions and take the smallest n number
131 // such that 2^n < m. n it will be order of the hilbert curve
132
133 size_t max = 0;
134
135 for (size_t i = 0; i < dim ; i++)
136 {
137 if (max < gr.size(i))
138 max = gr.size(i);
139 }
140
141 // Get the order of the hilbert-curve
142 size_t order = openfpm::math::log2_64(max);
143 if (1ul << order < max)
144 order += 1;
145
146 size_t n = 1 << order;
147
148 // Create the CellDecomoser
149
150 CellDecomposer_sm<dim,T> cd_sm;
151 cd_sm.setDimensions(domain, gr.getSize(), 0);
152
153 // create the hilbert curve
154
155 //hilbert curve iterator
156 grid_key_dx_iterator_hilbert<dim> h_it(order);
157
158 T spacing[dim];
159
160 // Calculate the hilbert curve spacing
161 for (size_t i = 0 ; i < dim ; i++)
162 spacing[i] = (domain.getHigh(i) - domain.getLow(i)) / n;
163
164 // Small grid to detect already assigned sub-sub-domains
165 grid_cpu<dim,aggregate<long int>> g(gr.getSize());
166 g.setMemory();
167
168 // Reset the grid to -1
169 grid_key_dx_iterator<dim> it(gr);
170 while (it.isNext())
171 {
172 auto key = it.get();
173
174 g.template get<0>(key) = -1;
175
176 ++it;
177 }
178
179 // Go along the hilbert-curve and divide the space
180
181 size_t proc_d = 0;
182 size_t ele_d = 0;
183 while (h_it.isNext())
184 {
185 auto key = h_it.get();
186
187 // Point p
188 Point<dim,T> p;
189
190 for (size_t i = 0 ; i < dim ; i++)
191 p.get(i) = key.get(i) * spacing[i] + spacing[i] / 2;
192
193 grid_key_dx<dim> sp = cd_sm.getCellGrid(p);
194
195 if (g.template get<0>(sp) == -1)
196 {
197 g.template get<0>(sp) = proc_d;
198 ele_d++;
199
200 if (ele_d >= accu.get(proc_d))
201 proc_d++;
202 }
203
204 ++h_it;
205 }
206
207 // Fill from the grid to the graph
208
209 // Reset the grid to -1
210 grid_key_dx_iterator<dim> it2(gr);
211 while (it2.isNext())
212 {
213 auto key = it2.get();
214
215 gp.template vertex_p<nm_v_proc_id>(gr.LinId(key)) = g.template get<0>(key);
216
217 ++it2;
218 }
219
220 return;
221 }
222
223 /*! \brief Refine current decomposition
224 *
225 * Has no effect in this case
226 *
227 */
228 void refine()
229 {
230 std::cout << __FILE__ << ":" << __LINE__ << " You are trying to dynamicaly balance a fixed decomposition, this operation has no effect" << std::endl;
231 }
232
233 /*! \brief Compute the unbalance of the processor compared to the optimal balance
234 *
235 * \return the unbalance from the optimal one 0.01 mean 1%
236 */
237 float getUnbalance()
238 {
239 return gr.size() % v_cl.getProcessingUnits();
240 }
241
242 /*! \brief function that return the position of the vertex in the space
243 *
244 * \param id vertex id
245 * \param pos vector that will contain x, y, z
246 *
247 */
248 void getSubSubDomainPosition(size_t id, T (&pos)[dim])
249 {
250#ifdef SE_CLASS1
251 if (id >= gp.getNVertex())
252 std::cerr << __FILE__ << ":" << __LINE__ << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex() << ")\n";
253#endif
254
255 // Copy the geometrical informations inside the pos vector
256 pos[0] = gp.vertex(id).template get<nm_v_x>()[0];
257 pos[1] = gp.vertex(id).template get<nm_v_x>()[1];
258 if (dim == 3)
259 pos[2] = gp.vertex(id).template get<nm_v_x>()[2];
260 }
261
262 /*! \brief Function that set the weight of the vertex
263 *
264 * \param id vertex id
265 * \param weight to give to the vertex
266 *
267 */
268 inline void setComputationCost(size_t id, size_t weight)
269 {
270 std::cout << __FILE__ << ":" << __LINE__ << " You are trying to set the computation cost on a fixed decomposition, this operation has no effect" << std::endl;
271 }
272
273 /*! \brief Checks if weights are used on the vertices
274 *
275 * \return true if weights are used in the decomposition
276 */
277 bool weightsAreUsed()
278 {
279 return false;
280 }
281
282 /*! \brief function that get the weight of the vertex
283 *
284 * \param id vertex id
285 *
286 * \return the weight of the vertex
287 *
288 */
289 size_t getSubSubDomainComputationCost(size_t id)
290 {
291 return 1.0;
292 }
293
294 /*! \brief Compute the processor load counting the total weights of its vertices
295 *
296 * \return the computational load of the processor graph
297 */
298 size_t getProcessorLoad()
299 {
300 // Get the number of processing units
301 size_t Np = v_cl.getProcessingUnits();
302
303 // Calculate the best number of sub-domains for each
304 // processor
305 size_t N_tot = gr.size();
306 size_t N_best_each = N_tot / Np;
307 size_t N_rest = N_tot % Np;
308
309 if (v_cl.getProcessUnitID() < N_rest)
310 N_best_each += 1;
311
312 return N_best_each;
313 }
314
315 /*! \brief Set migration cost of the vertex id
316 *
317 * \param id of the vertex to update
318 * \param migration cost of the migration
319 */
320 void setMigrationCost(size_t id, size_t migration)
321 {
322 }
323
324 /*! \brief Set communication cost of the edge id
325 *
326 * \param v_id Id of the source vertex of the edge
327 * \param e i child of the vertex
328 * \param communication Communication value
329 */
330 void setCommunicationCost(size_t v_id, size_t e, size_t communication)
331 {
332 }
333
334 /*! \brief Returns total number of sub-sub-domains in the distribution graph
335 *
336 * \return number of sub-sub-domain
337 *
338 */
339 size_t getNSubSubDomains()
340 {
341 return gp.getNVertex();
342 }
343
344 /*! \brief Returns total number of neighbors of the sub-sub-domain id
345 *
346 * \param id id of the sub-sub-domain
347 */
348 size_t getNSubSubDomainNeighbors(size_t id)
349 {
350 return gp.getNChilds(id);
351 }
352
353 /*! \brief Print the current distribution and save it to VTK file
354 *
355 * \param file filename
356 *
357 */
358 void write(const std::string & file)
359 {
360 VTKWriter<Graph_CSR<nm_v<dim>, nm_e>, VTK_GRAPH> gv2(gp);
361 gv2.write(std::to_string(v_cl.getProcessUnitID()) + "_" + file + ".vtk");
362 }
363
364 const SpaceDistribution<dim,T> & operator=(const SpaceDistribution<dim,T> & dist)
365 {
366 gr = dist.gr;
367 domain = dist.domain;
368 gp = dist.gp;
369
370 return *this;
371 }
372
373 const SpaceDistribution<dim,T> & operator=(SpaceDistribution<dim,T> && dist)
374 {
375 v_cl = dist.v_cl;
376 gr = dist.gr;
377 domain = dist.domain;
378 gp.swap(dist.gp);
379
380 return *this;
381 }
382
383 /*! \brief It return the decomposition id
384 *
385 * It just return 0
386 *
387 * \return 0
388 *
389 */
390 size_t get_ndec()
391 {
392 return 0;
393 }
394};
395
396
397#endif /* SRC_DECOMPOSITION_DISTRIBUTION_SPACEDISTRIBUTION_HPP_ */
398