SpaceDistribution.hpp 8.18 KB
Newer Older
incardon's avatar
incardon committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337
/*
 * SpaceDistribution.hpp
 *
 *  Created on: Dec 3, 2016
 *      Author: i-bird
 */

#ifndef SRC_DECOMPOSITION_DISTRIBUTION_SPACEDISTRIBUTION_HPP_
#define SRC_DECOMPOSITION_DISTRIBUTION_SPACEDISTRIBUTION_HPP_

#include "util/mathutil.hpp"
#include "NN/CellList/CellDecomposer.hpp"

/*! \brief Class that distribute sub-sub-domains across processors using an hilbert curve
 *         to divide the space
 *
 * ### Initialize a Cartesian graph and decompose
 * \snippet Distribution_unit_tests.hpp Initialize a Space Cartesian graph and decompose
 *
 *
 */
template<unsigned int dim, typename T>
class SpaceDistribution
{
	//! Vcluster
	Vcluster & v_cl;

	//! Structure that store the cartesian grid information
	grid_sm<dim, void> gr;

	//! rectangular domain to decompose
	Box<dim, T> domain;

	//! Global sub-sub-domain graph
	Graph_CSR<nm_v, nm_e> gp;


public:

	/*! Constructor
	 *
	 * \param v_cl Vcluster to use as communication object in this class
	 */
	SpaceDistribution(Vcluster & v_cl)
	:v_cl(v_cl)
	{
	}

	/*! Copy constructor
	 *
	 * \param pm Distribution to copy
	 *
	 */
	SpaceDistribution(const ParMetisDistribution<dim,T> & pm)
	:v_cl(pm.v_cl)
	{
		this->operator=(pm);
	}

	/*! Copy constructor
	 *
	 * \param pm Distribution to copy
	 *
	 */
	SpaceDistribution(SpaceDistribution<dim,T> && pm)
	{
		this->operator=(pm);
	}

	/*! \brief Create the Cartesian graph
	 *
	 * \param grid info
	 * \param dom domain
	 */
	void createCartGraph(grid_sm<dim, void> & grid, Box<dim, T> dom)
	{
		size_t bc[dim];

		for (size_t i = 0 ; i < dim ; i++)
			bc[i] = NON_PERIODIC;

		// Set grid and domain
		gr = grid;
		domain = dom;

		// Create a cartesian grid graph
		CartesianGraphFactory<dim, Graph_CSR<nm_v, nm_e>> g_factory_part;
		gp = g_factory_part.template construct<NO_EDGE, nm_v::id, T, dim - 1, 0>(gr.getSize(), domain, bc);

		// Init to 0.0 axis z (to fix in graphFactory)
		if (dim < 3)
		{
			for (size_t i = 0; i < gp.getNVertex(); i++)
				gp.vertex(i).template get<nm_v::x>()[2] = 0.0;
		}
		for (size_t i = 0; i < gp.getNVertex(); i++)
			gp.vertex(i).template get<nm_v::global_id>() = i;

	}

	/*! \brief Get the current graph (main)
	 *
	 */
	Graph_CSR<nm_v, nm_e> & getGraph()
	{
		return gp;
	}

	/*! \brief Create the decomposition
	 *
	 */
	void decompose()
	{
		// Get the number of processing units
		size_t Np = v_cl.getProcessingUnits();

		// Calculate the best number of sub-domains for each
		// processor
		size_t N_tot = gr.size();
		size_t N_best_each = N_tot / Np;
		size_t N_rest = N_tot % Np;

		openfpm::vector<size_t> accu(Np);
		accu.get(0) = N_best_each + ((0 < N_rest)?1:0);
		for (size_t i = 1 ; i < Np ; i++)
			accu.get(i) = accu.get(i-1) + N_best_each + ((i < N_rest)?1:0);

		// Get the maximum along dimensions and take the smallest n number
		// such that 2^n < m. n it will be order of the hilbert curve

		size_t max = 0;

		for (size_t i = 0; i < dim ; i++)
		{
			if (max < gr.size(i))
				max = gr.size(i);
		}

		// Get the order of the hilbert-curve
		size_t order = openfpm::math::log2_64(max);
		if (1ul << order < max)
			order += 1;

		size_t n = 1 << order;

		// Create the CellDecomoser

		CellDecomposer_sm<dim,T> cd_sm;
		cd_sm.setDimensions(domain, gr.getSize(), 0);

		// create the hilbert curve

		//hilbert curve iterator
		grid_key_dx_iterator_hilbert<dim> h_it(order);

		T spacing[dim];

		// Calculate the hilbert curve spacing
		for (size_t i = 0 ; i < dim ; i++)
			spacing[i] = (domain.getHigh(i) - domain.getLow(i)) / n;

		// Small grid to detect already assigned sub-sub-domains
		grid_cpu<dim,aggregate<long int>> g(gr.getSize());
		g.setMemory();

		// Reset the grid to -1
		grid_key_dx_iterator<dim> it(gr);
		while (it.isNext())
		{
			auto key = it.get();

			g.template get<0>(key) = -1;

			++it;
		}

		// Go along the hilbert-curve and divide the space

		size_t proc_d = 0;
		size_t ele_d = 0;
		while (h_it.isNext())
		{
		  auto key = h_it.get();

		  // Point p
		  Point<dim,T> p;

		  for (size_t i = 0 ; i < dim ; i++)
			  p.get(i) = key.get(i) * spacing[i] + spacing[i] / 2;

		  grid_key_dx<dim> sp = cd_sm.getCellGrid(p);

		  if (g.template get<0>(sp) == -1)
		  {
			  g.template get<0>(sp) = proc_d;
			  ele_d++;

			  if (ele_d >= accu.get(proc_d))
				  proc_d++;
		  }

		  ++h_it;
		}

		// Fill from the grid to the graph

		// Reset the grid to -1
		grid_key_dx_iterator<dim> it2(gr);
		while (it2.isNext())
		{
			auto key = it2.get();

			gp.template vertex_p<nm_v::proc_id>(gr.LinId(key)) = g.template get<0>(key);

			++it2;
		}

		return;
	}

	/*! \brief Refine current decomposition
	 *
	 * Has no effect in this case
	 *
	 */
	void refine()
	{
		std::cout << __FILE__ << ":" << __LINE__ << " You are trying to dynamicaly balance a fixed decomposition, this operation has no effect" << std::endl;
	}

	/*! \brief Compute the unbalance of the processor compared to the optimal balance
	 *
	 * \return the unbalance from the optimal one 0.01 mean 1%
	 */
	float getUnbalance()
	{
		return gr.size() % v_cl.getProcessingUnits();
	}

	/*! \brief function that return the position of the vertex in the space
	 *
	 * \param id vertex id
	 * \param pos vector that will contain x, y, z
	 *
	 */
	void getSubSubDomainPosition(size_t id, T (&pos)[dim])
	{
		if (id >= gp.getNVertex())
			std::cerr << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex() << ")\n";

		// Copy the geometrical informations inside the pos vector
		pos[0] = gp.vertex(id).template get<nm_v::x>()[0];
		pos[1] = gp.vertex(id).template get<nm_v::x>()[1];
		if (dim == 3)
			pos[2] = gp.vertex(id).template get<nm_v::x>()[2];
	}

	/*! \brief Function that set the weight of the vertex
	 *
	 * \param id vertex id
	 * \param weight to give to the vertex
	 *
	 */
	inline void setComputationCost(size_t id, size_t weight)
	{
		std::cout << __FILE__ << ":" << __LINE__ << " You are trying to set the computation cost on a fixed decomposition, this operation has no effect" << std::endl;
	}

	/*! \brief Checks if weights are used on the vertices
	 *
	 * \return true if weights are used in the decomposition
	 */
	bool weightsAreUsed()
	{
		return false;
	}

	/*! \brief function that get the weight of the vertex
	 *
	 * \param id vertex id
	 *
	 */
	size_t getSubSubDomainComputationCost(size_t id)
	{
		return 1.0;
	}

	/*! \brief Compute the processor load counting the total weights of its vertices
	 *
	 * \return the computational load of the processor graph
	 */
	size_t getProcessorLoad()
	{
		// Get the number of processing units
		size_t Np = v_cl.getProcessingUnits();

		// Calculate the best number of sub-domains for each
		// processor
		size_t N_tot = gr.size();
		size_t N_best_each = N_tot / Np;
		size_t N_rest = N_tot % Np;

		if (v_cl.getProcessUnitID() < N_rest)
			N_best_each += 1;

		return N_best_each;
	}

	/*! \brief Set migration cost of the vertex id
	 *
	 * \param id of the vertex to update
	 * \param migration cost of the migration
	 */
	void setMigrationCost(size_t id, size_t migration)
	{
	}

	/*! \brief Set communication cost of the edge id
	 *
	 * \param v_id Id of the source vertex of the edge
	 * \param e i child of the vertex
	 * \param communication Communication value
	 */
	void setCommunicationCost(size_t v_id, size_t e, size_t communication)
	{
	}

	/*! \brief Returns total number of sub-sub-domains in the distribution graph
	 *
	 */
	size_t getNSubSubDomains()
	{
		return gp.getNVertex();
	}

	/*! \brief Returns total number of neighbors of the sub-sub-domain id
	 *
incardon's avatar
incardon committed
338
	 * \param id id of the sub-sub-domain
incardon's avatar
incardon committed
339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377
	 */
	size_t getNSubSubDomainNeighbors(size_t id)
	{
		return gp.getNChilds(id);
	}

	/*! \brief Print the current distribution and save it to VTK file
	 *
	 * \param file filename
	 *
	 */
	void write(const std::string & file)
	{
		VTKWriter<Graph_CSR<nm_v, nm_e>, VTK_GRAPH> gv2(gp);
		gv2.write(std::to_string(v_cl.getProcessUnitID()) + "_" + file + ".vtk");
	}

	const SpaceDistribution<dim,T> & operator=(const SpaceDistribution<dim,T> & dist)
	{
		gr = dist.gr;
		domain = dist.domain;
		gp = dist.gp;

		return *this;
	}

	const SpaceDistribution<dim,T> & operator=(SpaceDistribution<dim,T> && dist)
	{
		v_cl = dist.v_cl;
		gr = dist.gr;
		domain = dist.domain;
		gp.swap(dist.gp);

		return *this;
	}
};


#endif /* SRC_DECOMPOSITION_DISTRIBUTION_SPACEDISTRIBUTION_HPP_ */