CartDecomposition with Parmetis and Metis distribution

src/Decomposition/BasicDecomposition_unit_test.hpp

-#ifndef BASICDECOMPOSITION_UNIT_TEST_HPP
-#define BASICDECOMPOSITION_UNIT_TEST_HPP
-
-#include "BasicDecomposition.hpp"
-#include "util/mathutil.hpp"
-#include <boost/algorithm/string.hpp>
-
-BOOST_AUTO_TEST_SUITE (BasicDecomposition_test)
-
-#define SUB_UNIT_FACTOR 128
-#define DIM 2
-
-BOOST_AUTO_TEST_CASE( BasicDecomposition_test_use)
-{
-	// Vcluster
-	Vcluster & vcl = *global_v_cluster;
-
-	// Initialize the global VCluster
-	init_global_v_cluster(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv);
-
-	//! [Create CartDecomposition]
-	BasicDecomposition<DIM, float> dec(vcl);
-
-	// Physical domain
-	Box<DIM, float> box( { 0.0, 0.0, 0.0 }, { 1.0, 1.0, 1.0 });
-	size_t div[DIM];
-
-	// Get the number of processor and calculate the number of sub-domain
-	// for each processor (SUB_UNIT_FACTOR=64)
-	size_t n_proc = vcl.getProcessingUnits();
-	size_t n_sub = n_proc * SUB_UNIT_FACTOR;
-
-	// Set the number of sub-domains on each dimension (in a scalable way)
-	for (int i = 0; i < DIM; i++)
-	{
-		div[i] = openfpm::math::round_big_2(pow(n_sub,1.0/DIM));
-	}
-
-	// Define ghost
-	Ghost<DIM, float> g(0.01);
-
-	// Decompose
-	dec.setParameters(div, box, g);
-
-	//Add weights
-	float c_x = 0, c_y = 0, c_z = 0, radius2, eq;
-	float x = 0, y = 0, z = 0;
-	float thickness = 2;
-	int weight = 5;
-
-	size_t n_v = pow(div[0], DIM);
-
-	// Set center and radius of the distribution shape (Sphere)
-	radius2 = pow(2,2);
-	c_x = 2;
-	c_y = 2;
-
-	if (DIM == 3)
-		c_z = 2;
-
-	openfpm::vector<real_t> pos(DIM);
-	for (int i = 0; i < n_v; i++)
-	{
-		dec.getVertexPosition(i, pos);
-		x = pos.get(0) * 10;
-		y = pos.get(1) * 10;
-
-		if (DIM == 3)
-		z = pos.get(2) * 10;
-
-		eq = pow((x - c_x), 2) + pow((y - c_y), 2) + pow((z - c_z), 2);
-
-		if (eq <= radius2)
-		{
-			dec.setVertexWeight(i, weight);
-		} else
-		{
-			dec.setVertexWeight(i, 1);
-		}
-	}
-
-	dec.decompose();
-
-	dec.printCurrentGraph(0);
-
-	float stime = 0.1, etime = 5, tstep = 0.1;
-
-	for(real_t t = stime, i = 1, t_sim = 1;
-			t < etime;
-			t = t + tstep, i++, t_sim++)
-	{
-		clock_t start = clock();
-
-		if(t < etime/2){
-			c_x += tstep;
-			c_y += tstep;
-		}else{
-			c_x -= tstep;
-			c_y -= tstep;
-		}
-
-
-		if (DIM == 3)
-			c_z += tstep;
-
-		openfpm::vector<real_t> pos(DIM);
-		for (int i = 0; i < n_v; i++)
-		{
-			dec.getVertexPosition(i, pos);
-			x = pos.get(0) * 10;
-			y = pos.get(1) * 10;
-
-			if (DIM == 3)
-			z = pos.get(2) * 10;
-
-			eq = pow((x - c_x), 2) + pow((y - c_y), 2) + pow((z - c_z), 2);
-
-			if (eq <= radius2)
-			{
-				dec.setVertexWeight(i, weight);
-			} else
-			{
-				dec.setVertexWeight(i, 1);
-			}
-		}
-
-		usleep(1000*t_sim);
-
-		clock_t end = clock();
-		long time = end - start;
-
-		bool needed = dec.balanceNeeded(time);
-
-		if(needed){
-			dec.refine();
-			t_sim = 1;
-		}
-
-		dec.printCurrentGraph(i);
-	}
-
-	//print statistics
-	if(vcl.getProcessUnitID() == 0)
-	{
-		float avg = dec.getTotalMovedV()/((etime-stime)/tstep);
-
-		std::cout << "Moved vertices average: " << avg << "\n";
-		std::cout << "Max number of moved vertices: " << dec.getMaxMovedV() << "\n";
-	}
-
-	// create a ghost border
-	dec.calculateGhostBoxes();
-
-	// For each calculated ghost box
-	for (size_t i = 0; i < dec.getNIGhostBox(); i++)
-	{
-		SpaceBox<DIM,float> b = dec.getIGhostBox(i);
-		size_t proc = dec.getIGhostBoxProcessor(i);
-
-		// sample one point inside the box
-		Point<DIM,float> p = b.rnd();
-
-		// Check that ghost_processorsID return that processor number
-		const openfpm::vector<size_t> & pr = dec.template ghost_processorID<BasicDecomposition<DIM,float>::processor_id>(p);
-
-		bool found = false;
-
-		for (size_t j = 0; j < pr.size(); j++)
-		{
-			if (pr.get(j) == proc)
-			{	found = true; break;}
-		}
-
-		if (found == false)
-		{
-			const openfpm::vector<size_t> pr2 = dec.template ghost_processorID<BasicDecomposition<DIM,float>::processor_id>(p);
-		}
-
-		BOOST_REQUIRE_EQUAL(found,true);
-	}
-
-	// Check the consistency
-
-	bool val = dec.check_consistency();
-	BOOST_REQUIRE_EQUAL(val,true);
-}
-
-BOOST_AUTO_TEST_SUITE_END()
-
-#endif

src/Decomposition/CartDecomposition_unit_test.hpp

@@ -3,10 +3,13 @@
 
 #include "CartDecomposition.hpp"
 #include "util/mathutil.hpp"
+#include "DLB.hpp"
+#include <boost/algorithm/string.hpp>
 
-BOOST_AUTO_TEST_SUITE( CartDecomposition_test )
+BOOST_AUTO_TEST_SUITE (CartDecomposition_test)
 
-#define SUB_UNIT_FACTOR 64
+#define SUB_UNIT_FACTOR 128
+#define DIM 2
 
 BOOST_AUTO_TEST_CASE( CartDecomposition_test_use)
 {
@@ -17,11 +20,14 @@ BOOST_AUTO_TEST_CASE( CartDecomposition_test_use)
 	init_global_v_cluster(&boost::unit_test::framework::master_test_suite().argc,&boost::unit_test::framework::master_test_suite().argv);
 
 	//! [Create CartDecomposition]
-	CartDecomposition<3,float> dec(vcl);
+	CartDecomposition<DIM, float> dec(vcl);
+
+	// Init DLB tool
+	DLB dlb(vcl);
 
 	// Physical domain
-	Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
-	size_t div[3];
+	Box<DIM, float> box( { 0.0, 0.0, 0.0 }, { 1.0, 1.0, 1.0 });
+	size_t div[DIM];
 
 	// Get the number of processor and calculate the number of sub-domain
 	// for each processor (SUB_UNIT_FACTOR=64)
@@ -29,43 +35,157 @@ BOOST_AUTO_TEST_CASE( CartDecomposition_test_use)
 	size_t n_sub = n_proc * SUB_UNIT_FACTOR;
 
 	// Set the number of sub-domains on each dimension (in a scalable way)
-	for (int i = 0 ; i < 3 ; i++)
-	{div[i] = openfpm::math::round_big_2(pow(n_sub,1.0/3));}
+	for (int i = 0; i < DIM; i++)
+	{
+		div[i] = openfpm::math::round_big_2(pow(n_sub,1.0/DIM));
+	}
 
 	// Define ghost
-	Ghost<3,float> g(0.01);
+	Ghost<DIM, float> g(0.01);
 
 	// Decompose
-	dec.setParameters(div,box,g);
+	dec.setParameters(div, box, g);
+
+	//Add weights
+	float c_x = 0, c_y = 0, c_z = 0, radius2, eq;
+	float x = 0, y = 0, z = 0;
+	float thickness = 2;
+	int weight = 5;
+
+	size_t n_v = pow(div[0], DIM);
+
+	// Set center and radius of the distribution shape (Sphere)
+	radius2 = pow(2,2);
+	c_x = 2;
+	c_y = 2;
+
+	if (DIM == 3)
+	c_z = 2;
+
+	openfpm::vector<real_t> pos(DIM);
+	for (int i = 0; i < n_v; i++)
+	{
+		dec.getSubSubDomainPosition(i, pos);
+		x = pos.get(0) * 10;
+		y = pos.get(1) * 10;
+
+		if (DIM == 3)
+		z = pos.get(2) * 10;
+
+		eq = pow((x - c_x), 2) + pow((y - c_y), 2) + pow((z - c_z), 2);
+
+		if (eq <= radius2)
+		{
+			dec.setSubSubDomainComputationCost(i, weight);
+		}
+		else
+		{
+			dec.setSubSubDomainComputationCost(i, 1);
+		}
+	}
+
+	dec.decompose();
+
+	dec.printCurrentDecomposition(0);
+
+	float stime = 0.1, etime = 5, tstep = 0.1;
+
+	dlb.setSimulationStartTime(0);
+	dlb.setSimulationEndTime(5);
+
+	for(real_t t = stime, i = 1, t_sim = 1;
+			t < etime;
+			t = t + tstep, i++, t_sim++)
+	{
+		dlb.setIterationStartTime(clock());
+
+		if(t < etime/2)
+		{
+			c_x += tstep;
+			c_y += tstep;
+		}
+		else
+		{
+			c_x -= tstep;
+			c_y -= tstep;
+		}
+
+		if (DIM == 3)
+		c_z += tstep;
+
+		openfpm::vector<real_t> pos(DIM);
+		for (int i = 0; i < n_v; i++)
+		{
+			dec.getSubSubDomainPosition(i, pos);
+			x = pos.get(0) * 10;
+			y = pos.get(1) * 10;
+
+			if (DIM == 3)
+			z = pos.get(2) * 10;
+
+			eq = pow((x - c_x), 2) + pow((y - c_y), 2) + pow((z - c_z), 2);
+
+			if (eq <= radius2)
+			{
+				dec.setSubSubDomainComputationCost(i, weight);
+			}
+			else
+			{
+				dec.setSubSubDomainComputationCost(i, 1);
+			}
+		}
+
+		usleep(1000*t_sim);
+
+		dlb.setIterationEndTime(clock());
+
+		dec.rebalance(dlb);
+
+		if(dlb.rebalanceNeeded())
+		{
+			t_sim = 1;
+		}
+
+		dec.printCurrentDecomposition(i);
+	}
+
+	//print statistics
+	/*
+	 if(vcl.getProcessUnitID() == 0)
+	 {
+	 float avg = dec.getTotalMovedV()/((etime-stime)/tstep);
+
+	 std::cout << "Moved vertices average: " << avg << "\n";
+	 std::cout << "Max number of moved vertices: " << dec.getMaxMovedV() << "\n";
+	 }
+	 */
 
 	// create a ghost border
 	dec.calculateGhostBoxes();
 
-	//! [Create CartDecomposition]
-
 	// For each calculated ghost box
-	for (size_t i = 0 ; i < dec.getNIGhostBox() ; i++)
+	for (size_t i = 0; i < dec.getNIGhostBox(); i++)
 	{
-		SpaceBox<3,float> b = dec.getIGhostBox(i);
+		SpaceBox<DIM,float> b = dec.getIGhostBox(i);
 		size_t proc = dec.getIGhostBoxProcessor(i);
 
 		// sample one point inside the box
-		Point<3,float> p = b.rnd();
+		Point<DIM,float> p = b.rnd();
 
 		// Check that ghost_processorsID return that processor number
-		const openfpm::vector<size_t> & pr = dec.template ghost_processorID<CartDecomposition<3,float>::processor_id>(p);
+		const openfpm::vector<size_t> & pr = dec.template ghost_processorID<CartDecomposition<DIM,float>::processor_id>(p);
 
 		bool found = false;
 
-		for (size_t j = 0; j < pr.size() ; j++)
+		for (size_t j = 0; j < pr.size(); j++)
 		{
 			if (pr.get(j) == proc)
-			{found = true; break;}
+			{	found = true; break;}
 		}
 
 		if (found == false)
 		{
-			const openfpm::vector<size_t> pr2 = dec.template ghost_processorID<CartDecomposition<3,float>::processor_id>(p);
+			const openfpm::vector<size_t> pr2 = dec.template ghost_processorID<CartDecomposition<DIM,float>::processor_id>(p);
 		}
 
 		BOOST_REQUIRE_EQUAL(found,true);
@@ -79,5 +199,4 @@ BOOST_AUTO_TEST_CASE( CartDecomposition_test_use)
 
 BOOST_AUTO_TEST_SUITE_END()
 
-
 #endif

src/Decomposition/MetisDistribution.hpp

@@ -2,14 +2,207 @@
  * MetisDistribution.hpp
  *
  *  Created on: Nov 19, 2015
- *      Author: i-bird
+ *      Author: Antonio Leo
  */
 
 #ifndef SRC_DECOMPOSITION_METISDISTRIBUTION_HPP_
 #define SRC_DECOMPOSITION_METISDISTRIBUTION_HPP_
 
+#include "SubdomainGraphNodes.hpp"
+#include "metis_util.hpp"
 
+template<unsigned int dim, typename T, template<unsigned int, typename > class Domain = Box>
+class MetisDistribution
+{
+	//! Vcluster
+	Vcluster & v_cl;
 
+	//! Structure that store the cartesian grid information
+	grid_sm<dim, void> gr;
 
+	//! rectangular domain to decompose
+	Domain<dim, T> domain;
+
+	//! Global sub-sub-domain graph
+	Graph_CSR<nm_v, nm_e> gp;
+
+	//! Flag to check if weights are used on vertices
+	bool verticesGotWeights = false;
+
+public:
+
+	//! constructor
+	MetisDistribution(Vcluster & v_cl) :
+			v_cl(v_cl)
+	{
+	}
+
+	/*! \brief Initialize the distribution graph
+	 *
+	 */
+	void init(grid_sm<dim, void> & grid, Domain<dim, T> dom)
+	{
+		// 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, 1, 2>(gr.getSize(), domain);
+
+		// 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::z>() = 0.0;
+			}
+		}
+	}
+
+	/*! \brief Get the current graph (main)
+	 *
+	 */
+	Graph_CSR<nm_v, nm_e> & getGraph()
+	{
+		return gp;
+	}
+
+	/*! \brief Create first decomposition, it divides the graph in slices and give each slice to a processor
+	 *
+	 */
+	void decompose()
+	{
+		Metis<Graph_CSR<nm_v, nm_e>> met(gp, v_cl.getProcessingUnits(), verticesGotWeights);
+
+		// decompose
+		met.decompose<nm_v::proc_id>();
+	}
+
+	/*! \brief Refine current decomposition (NOT AVAILABLE on Metis)
+	 *
+	 * It has no function
+	 */
+	void refine()
+	{
+	}
+
+	/*! \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 getVertexPosition(size_t id, openfpm::vector<real_t> &pos)
+	{
+		if (id >= gp.getNVertex())
+			std::cerr << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex()
+					<< ")\n";
+
+		pos.get(0) = gp.vertex(id).template get<nm_v::x>();
+		pos.get(1) = gp.vertex(id).template get<nm_v::y>();
+
+		if (dim == 3)
+			pos.get(2) = gp.vertex(id).template get<nm_v::z>();
+	}
+
+	/*! \brief function that set the weight of the vertex
+	 *
+	 * \param id vertex id
+	 * \param wieght to give to the vertex
+	 *
+	 */
+	void setVertexWeight(size_t id, size_t weight)
+	{
+		if(!verticesGotWeights)
+			verticesGotWeights = true;
+
+		if (id >= gp.getNVertex())
+			std::cerr << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex()
+					<< ")\n";
+
+		gp.vertex(id).template get<nm_v::computation>() = weight;
+	}
+
+	/*! \brief Checks if weights are used on the vertices
+	 *
+	 */
+	bool weightsAreUsed()
+	{
+		return verticesGotWeights;
+	}
+
+	/*! \brief function that get the weight of the vertex
+	 *
+	 * \param id vertex id
+	 *
+	 */
+	size_t getVertexWeight(size_t id)
+	{
+		if (id >= gp.getNVertex())
+			std::cerr << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex()
+					<< ")\n";
+
+		return gp.vertex(id).template get<nm_v::computation>();
+	}
+
+	/*! \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)
+	{
+		if (id >= gp.getNVertex())
+			std::cerr << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex()
+					<< ")\n";
+
+		gp.vertex(id).template get<nm_v::migration>() = migration;
+	}
+
+	/*! \brief Set communication cost of the edge id
+	 *
+	 */
+	void setCommunicationCost(size_t id, size_t communication)
+	{
+		if (id >= gp.getNEdge())
+			std::cerr << "Such edge doesn't exist (id = " << id << ", " << "total size = " << gp.getNEdge() << ")\n";
+
+		gp.edge(id).template get<nm_e::communication>() = 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
+	 *
+	 * \param i id of the sub-sub-domain
+	 */
+	size_t getNSubSubDomainNeighbors(size_t id)
+	{
+		if (id >= gp.getNVertex())
+			std::cerr << "Such vertex doesn't exist (id = " << id << ", " << "total size = " << gp.getNVertex()
+					<< ")\n";
+
+		return gp.getNChilds(id);
+	}
+
+	/* \brief Print current graph and save it to file with name test_graph_[id]
+	 *
+	 * \param id to attach to the filename
+	 *
+	 */
+	void printCurrentDecomposition(int id)
+	{
+		VTKWriter<Graph_CSR<nm_v, nm_e>, GRAPH> gv2(gp);
+		gv2.write("test_graph_" + std::to_string(id) + ".vtk");
+
+	}
+};
 
 #endif /* SRC_DECOMPOSITION_METISDISTRIBUTION_HPP_ */

src/main.cpp

@@ -25,5 +25,4 @@
 #include "dec_optimizer_unit_test.hpp"
 #include "Grid/grid_dist_id_unit_test.hpp"
 #include "Vector/vector_dist_unit_test.hpp"
-#include "Decomposition/BasicDecomposition_unit_test.hpp"
-
+#include "Decomposition/DLB_unit_test.hpp"

src/metis_util.hpp

@@ -11,6 +11,7 @@
 #include <iostream>
 #include "metis.h"
 #include "VTKWriter.hpp"