Adding interface for Multiphase DLB

example/Vector/0_simple/main.cpp

@@ -146,7 +146,7 @@ int main(int argc, char* argv[])
 
 	//! \cond [vector instantiation] \endcond
 
-	vector_dist<2,float, aggregate<float,float[3],float[3][3]> > vd(4096,domain,bc,g);
+        vector_dist<2,float, aggregate<float,float[3],float[3][3]> > vd(4096,domain,bc,g);
 
 	// the scalar is the element at position 0 in the aggregate
 	const int scalar = 0;
@@ -287,7 +287,7 @@ int main(int argc, char* argv[])
 		vd.template getProp<tensor>(p)[1][2] = 1.0;
 		vd.template getProp<tensor>(p)[2][0] = 1.0;
 		vd.template getProp<tensor>(p)[2][1] = 1.0;
-		vd.template getProp<tensor>(p)[2][2] = 1.0;
+        vd.template getProp<tensor>(p)[2][2] = 1.0;
 
 		// increment the counter
 		cnt++;

src/Decomposition/Distribution/SpaceDistribution.hpp

@@ -281,6 +281,8 @@ public:
 	 *
 	 * \param id vertex id
 	 *
+	 * \return the weight of the vertex
+	 *
 	 */
 	size_t getSubSubDomainComputationCost(size_t id)
 	{
@@ -328,6 +330,8 @@ public:
 	}
 
 	/*! \brief Returns total number of sub-sub-domains in the distribution graph
+	 *
+	 * \return number of sub-sub-domain
 	 *
 	 */
 	size_t getNSubSubDomains()

src/Grid/grid_dist_id.hpp

@@ -939,9 +939,9 @@ public:
 		return true;
 	}
 
-	/*! \brief Get the size of local domain grids
+	/*! \brief Get the total number of grid points for the calling processor
 	 *
-	 * \return The size of the local domain
+	 * \return The number of grid points
 	 *
 	 */
 	size_t getLocalDomainSize() const
@@ -959,9 +959,9 @@ public:
 		return total;
 	}
 
-	/*! \brief Get the size of local domain grids
+	/*! \brief Get the total number of grid points with ghost for the calling processor
 	 *
-	 * \return The size of the local domain
+	 * \return The number of grid points
 	 *
 	 */
 	size_t getLocalDomainWithGhostSize() const

src/Vector/se_class3_vector.hpp

@@ -38,17 +38,19 @@ enum ptype
 	INSIDE
 };
 
-// is initialized
+//! is initialized
 template<typename T>
 struct is_initialized
 {
+	//! it indicate the property is not initialized
 	static const int init = UNINITIALIZED;
 };
 
-// is initialized
+//! is initialized
 template<typename T>
 struct is_initialized<openfpm::vector<T>>
 {
+	//! it indicaste that property is clean
 	static const int init = CLEAN;
 };
 
@@ -61,22 +63,21 @@ struct is_initialized<openfpm::vector<T>>
  * element of the boost::vector the operator() is called.
  * Is mainly used to initialize the properties
  *
- * \tparam encap source
- * \tparam encap dst
+ * \tparam Np number of properties
+ * \tparam vector type of vector
  *
  */
 
 template<unsigned int Np, typename vector>
 struct init_prop
 {
-	//! vector for prop initializetion
+	//! vector for prop initialization
 	size_t (& prp_init)[Np];
 
 	/*! \brief constructor
 	 *
 	 *
-	 * \param src encapsulated object1
-	 * \param dst encapsulated object2
+	 * \param prp_init properties to initialize
 	 *
 	 */
 	inline init_prop(size_t ( & prp_init)[Np])
@@ -103,11 +104,13 @@ struct init_prop
 template<typename tcheck, bool foundamental>
 struct typeCheck
 {
+	//! It check if the type is Nan, data type to check
 	static bool isNan(const tcheck & data)
 	{
 		return false;
 	}
 
+	//! It check is the type is infinity, data type to checl
 	static bool isInf(const tcheck & data)
 	{
 		return false;
@@ -234,7 +237,7 @@ struct propCheckNAN
 	{
 		typedef typename boost::mpl::at<typename vector::value_type::type,typename boost::mpl::int_<T::value> >::type type_to_check;
 
-		bool snn = typeCheck<type_to_check,std::is_fundamental<type_to_check>::value>::isNan(data.template getProp<T::value>(id));
+		bool snn = typeCheck<type_to_check,std::is_fundamental<type_to_check>::value>::isNan(data.template getPropNC<T::value>(id));
 
 		if (snn == true)
 		{
@@ -266,7 +269,8 @@ struct propCheckINF
 
 	/*! \brief constructor
 	 *
-	 * \param
+	 * \param check the the property is infinity
+	 * \param id element
 	 *
 	 */
 	inline propCheckINF(const vector & data, size_t id)
@@ -311,6 +315,54 @@ static inline std::string getParticleTypeString(size_t type)
 	return std::string();
 }
 
+template<unsigned int prp, unsigned int Np, typename vector>  void check_for_pos_nan_inf(const vector & vd, size_t p)
+{
+#ifdef CHECKFOR_POSINF
+
+			if ( std::isinf(vd.getPosNC(p)[0]) || std::isinf(vd.getPosNC(p)[1]) || std::isinf(vd.getPosNC(p)[2]) )
+			{
+				std::cerr << __FILE__ << ":" << __LINE__ << " error detected INF in position for particle p=" << p << " of type=" << getParticleTypeString(vd.template getPropNC<Np+SE3_TYPE>(p)) << std::endl;
+				ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
+			}
+
+#endif
+
+#ifdef CHECKFOR_POSNAN
+
+			if ( std::isnan(vd.getPosNC(p)[0]) || std::isnan(vd.getPosNC(p)[1]) || std::isnan(vd.getPosNC(p)[2]) )
+			{
+				std::cerr << __FILE__ << ":" << __LINE__ << " error detected NAN in position for particle p=" << p << " of type=" << getParticleTypeString(vd.template getPropNC<Np+SE3_TYPE>(p)) << std::endl;
+				ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
+			}
+
+#endif
+}
+
+template<unsigned int prp, unsigned int Np_real, typename vector>  void check_for_prop_nan_inf(const vector & vd, size_t p)
+{
+#ifdef CHECKFOR_PROPINF
+
+			{
+				propCheckINF<vector> checker(vd,p);
+
+				boost::mpl::for_each_ref< boost::mpl::range_c<int,0, Np_real > > (checker);
+			}
+
+#endif
+
+#ifdef CHECKFOR_PROPNAN
+
+			{
+				propCheckNAN<vector> checker(vd,p);
+
+				boost::mpl::for_each_ref< boost::mpl::range_c<int,0, Np_real > >(checker);
+			}
+
+#endif
+}
+
+
+
 /*! \brief This class check for inconsistency access
  *
  * \tparam Np number of properties
@@ -707,46 +759,8 @@ class se_class3_vector
 				}
 			}
 
-#ifdef CHECK_FOR_POSINF
-
-			if ( std::isinf(vd.getPosNC(p)[0]) || std::isinf(vd.getPosNC(p)[1]) || std::isinf(vd.getPosNC(p)[2]) )
-			{
-				std::cerr << __FILE__ << ":" << __LINE__ << " error detected INF in position for particle p=" << p << " of type=" << getParticleTypeString(vd.template getPropNC<Np+SE3_TYPE>(p)) << std::endl;
-				ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
-			}
-
-#endif
-
-#ifdef CHECKFOR_POSNAN
-
-			if ( std::isnan(vd.getPosNC(p)[0]) || std::isnan(vd.getPosNC(p)[1]) || std::isnan(vd.getPosNC(p)[2]) )
-			{
-				std::cerr << __FILE__ << ":" << __LINE__ << " error detected NAN in position for particle p=" << p << " of type=" << getParticleTypeString(vd.template getPropNC<Np+SE3_TYPE>(p)) << std::endl;
-				ACTION_ON_ERROR(VECTOR_DIST_ERROR_OBJECT);
-			}
-
-#endif
-
-#ifdef CHECKFOR_PROPINF
-
-			{
-				propCheckINF<vector> checker(vd,p);
-
-				boost::mpl::for_each_ref< boost::mpl::range_c<int,0, Np_real > > (checker);
-			}
-
-#endif
-
-#ifdef CHECKFOR_PROPNAN
-
-			{
-				propCheckNAN<vector> checker(vd,p);
-
-				boost::mpl::for_each_ref< boost::mpl::range_c<int,0, Np_real > >(checker);
-			}
-
-#endif
-
+			check_for_pos_nan_inf<prp>(vd,p);
+			check_for_prop_nan_inf<prp,Np_real>(vd,p);
 		}
 
 		template<unsigned int prp> void write(vector & vd, size_t p)

src/Vector/vector_dist.hpp

@@ -231,6 +231,7 @@ private:
 		}
 	}
 
+
 public:
 
 	//! space type
@@ -424,6 +425,10 @@ public:
 	 */
 	inline auto getPos(vect_dist_key_dx vec_key) -> decltype(v_pos.template get<0>(vec_key.getKey()))
 	{
+#ifdef SE_CLASS3
+		check_for_pos_nan_inf<prop::max_prop_real,prop::max_prop>(*this,vec_key.getKey());
+#endif
+
 		return v_pos.template get<0>(vec_key.getKey());
 	}
 
@@ -438,6 +443,9 @@ public:
 	 */
 	inline auto getPos(vect_dist_key_dx vec_key) const -> decltype(v_pos.template get<0>(vec_key.getKey()))
 	{
+#ifdef SE_CLASS3
+		check_for_pos_nan_inf<prop::max_prop_real,prop::max_prop>(*this,vec_key.getKey());
+#endif
 		return v_pos.template get<0>(vec_key.getKey());
 	}
 
@@ -452,6 +460,9 @@ public:
 	 */
 	inline auto getPos(size_t vec_key) -> decltype(v_pos.template get<0>(vec_key))
 	{
+#ifdef SE_CLASS3
+		check_for_pos_nan_inf<prop::max_prop_real,prop::max_prop>(*this,vec_key);
+#endif
 		return v_pos.template get<0>(vec_key);
 	}
 
@@ -466,6 +477,9 @@ public:
 	 */
 	inline auto getPos(size_t vec_key) const -> decltype(v_pos.template get<0>(vec_key))
 	{
+#ifdef SE_CLASS3
+		check_for_pos_nan_inf<prop::max_prop_real,prop::max_prop>(*this,vec_key);
+#endif
 		return v_pos.template get<0>(vec_key);
 	}
 
@@ -481,6 +495,9 @@ public:
 	 */
 	template<unsigned int id> inline auto getProp(vect_dist_key_dx vec_key) -> decltype(v_prp.template get<id>(vec_key.getKey()))
 	{
+#ifdef SE_CLASS3
+		check_for_prop_nan_inf<id,prop::max_prop+SE3_STATUS>(*this,vec_key.getKey());
+#endif
 		return v_prp.template get<id>(vec_key.getKey());
 	}
 
@@ -496,6 +513,9 @@ public:
 	 */
 	template<unsigned int id> inline auto getProp(vect_dist_key_dx vec_key) const -> decltype(v_prp.template get<id>(vec_key.getKey()))
 	{
+#ifdef SE_CLASS3
+		check_for_prop_nan_inf<id,prop::max_prop+SE3_STATUS>(*this,vec_key.getKey());
+#endif
 		return v_prp.template get<id>(vec_key.getKey());
 	}
 
@@ -511,6 +531,9 @@ public:
 	 */
 	template<unsigned int id> inline auto getProp(size_t vec_key) -> decltype(v_prp.template get<id>(vec_key))
 	{
+#ifdef SE_CLASS3
+		check_for_prop_nan_inf<id,prop::max_prop+SE3_STATUS>(*this,vec_key);
+#endif
 		return v_prp.template get<id>(vec_key);
 	}
 
@@ -526,6 +549,9 @@ public:
 	 */
 	template<unsigned int id> inline auto getProp(size_t vec_key) const -> decltype(v_prp.template get<id>(vec_key))
 	{
+#ifdef SE_CLASS3
+		check_for_prop_nan_inf<id,prop::max_prop+SE3_STATUS>(*this,vec_key);
+#endif
 		return v_prp.template get<id>(vec_key);
 	}
 
@@ -1715,14 +1741,76 @@ public:
 	 * from the particles
 	 *
 	 * \param md Model to use
+	 * \param vd vector to add for the computational cost
 	 * \param ts It is an optional parameter approximately should be the number of ghost get between two
 	 *           rebalancing at first decomposition this number can be ignored (default = 1) because not used
 	 *
 	 */
-	template <typename Model=ModelLin>inline void addComputationCosts(Model md=Model(), size_t ts = 1)
+	template <typename Model=ModelLin>inline void addComputationCosts(const self & vd, Model md=Model())
 	{
 		CellDecomposer_sm<dim, St, shift<dim,St>> cdsm;
 
+		Decomposition & dec = getDecomposition();
+
+		cdsm.setDimensions(dec.getDomain(), dec.getDistGrid().getSize(), 0);
+
+		auto it = vd.getDomainIterator();
+
+		while (it.isNext())
+		{
+			size_t v = cdsm.getCell(vd.getPos(it.get()));
+
+			md.addComputation(dec,vd,v,it.get().getKey());
+
+			++it;
+		}
+	}
+
+	template <typename Model=ModelLin> void finalizeComputationCosts(Model md=Model(), size_t ts = 1)
+	{
+		Decomposition & dec = getDecomposition();
+		auto & dist = getDecomposition().getDistribution();
+
+		dec.computeCommunicationAndMigrationCosts(ts);
+
+		// Go throught all the sub-sub-domains and apply the model
+
+		for (size_t i = 0 ; i < dist.getNOwnerSubSubDomains(); i++)
+		{md.applyModel(dec,dist.getOwnerSubSubDomain(i));}
+
+		dist.setDistTol(md.distributionTol());
+	}
+
+	/*! \brief Initialize the computational cost
+	 *
+	 */
+	void initializeComputationCosts()
+	{
+		Decomposition & dec = getDecomposition();
+		auto & dist = getDecomposition().getDistribution();
+
+		for (size_t i = 0; i < dist.getNOwnerSubSubDomains() ; i++)
+		{dec.setSubSubDomainComputationCost(dist.getOwnerSubSubDomain(i) , 1);}
+	}
+
+	/*! \brief Add the computation cost on the decomposition coming
+	 * from the particles
+	 *
+	 * \param md Model to use
+	 * \param ts It is an optional parameter approximately should be the number of ghost get between two
+	 *           rebalancing at first decomposition this number can be ignored (default = 1) because not used
+	 *
+	 */
+	template <typename Model=ModelLin>inline void addComputationCosts(Model md=Model(), size_t ts = 1)
+	{
+		initializeComputationCosts();
+
+		addComputationCosts(*this,md);
+
+		finalizeComputationCosts(md,ts);
+
+/*		CellDecomposer_sm<dim, St, shift<dim,St>> cdsm;
+
 		Decomposition & dec = getDecomposition();
 		auto & dist = getDecomposition().getDistribution();
 
@@ -1749,7 +1837,7 @@ public:
 		for (size_t i = 0 ; i < dist.getNOwnerSubSubDomains(); i++)
 			md.applyModel(dec,dist.getOwnerSubSubDomain(i));
 
-		dist.setDistTol(md.distributionTol());
+		dist.setDistTol(md.distributionTol());*/
 	}
 
 	/*! \brief Save the distributed vector on HDF5 file
@@ -1996,7 +2084,7 @@ public:
 		sz = 0;
 
 		for (size_t i = 0 ; i < v_cl.getProcessUnitID() ; i++)
-			sz += accu.get(i);
+		{sz += accu.get(i);}
 
 		return sz;
 	}

src/Vector/vector_dist_dlb_test.hpp

@@ -46,8 +46,6 @@ template<typename vector_type> void test_dlb_vector()
 	vd.getDecomposition().decompose();
 	vd.map();
 
-	vd.getDecomposition().getDistribution().write("dist_out_");
-	vd.getDecomposition().write("dec_out_");
 
 	vd.addComputationCosts(md);
 
@@ -112,11 +110,226 @@ template<typename vector_type> void test_dlb_vector()
 	}
 }
 
+
+template<typename vector_type> void test_dlb_multi_phase_vector()
+{
+	Vcluster & v_cl = create_vcluster();
+
+	if (v_cl.getProcessingUnits() > 8)
+		return;
+
+	Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
+	Ghost<3,float> g(0.1);
+	size_t bc[3] = {PERIODIC,PERIODIC,PERIODIC};
+
+	vector_type vd0(0,domain,bc,g,DEC_GRAN(2048));
+	vector_type vd1(0,domain,bc,g,DEC_GRAN(2048));
+	vector_type vd2(0,domain,bc,g,DEC_GRAN(2048));
+	vector_type vd3(0,domain,bc,g,DEC_GRAN(2048));
+
+	// Only processor 0 initialy add particles on a corner of a domain
+
+	if (v_cl.getProcessUnitID() == 0)
+	{
+		for(size_t i = 0 ; i < 50000 ; i++)
+		{
+			vd0.add();
+			vd1.add();
+			vd2.add();
+			vd3.add();
+
+			vd0.getLastPos()[0] = ((float)rand())/RAND_MAX * 0.3;
+			vd0.getLastPos()[1] = ((float)rand())/RAND_MAX * 0.3;
+			vd0.getLastPos()[2] = ((float)rand())/RAND_MAX * 0.3;
+
+			vd1.getLastPos()[0] = ((float)rand())/RAND_MAX * 0.3 + 0.1;
+			vd1.getLastPos()[1] = ((float)rand())/RAND_MAX * 0.3 + 0.1;
+			vd1.getLastPos()[2] = ((float)rand())/RAND_MAX * 0.3 + 0.1;
+
+			vd2.getLastPos()[0] = ((float)rand())/RAND_MAX * 0.3 + 0.2;
+			vd2.getLastPos()[1] = ((float)rand())/RAND_MAX * 0.3 + 0.2;
+			vd2.getLastPos()[2] = ((float)rand())/RAND_MAX * 0.3 + 0.2;
+
+			vd3.getLastPos()[0] = ((float)rand())/RAND_MAX * 0.3 + 0.3;
+			vd3.getLastPos()[1] = ((float)rand())/RAND_MAX * 0.3 + 0.3;
+			vd3.getLastPos()[2] = ((float)rand())/RAND_MAX * 0.3 + 0.3;
+		}
+	}
+
+	vd0.map();
+	vd0.template ghost_get<>();
+	vd1.map();
+	vd1.template ghost_get<