Before refactor communication of vector dist

CHANGELOG.md

 # Change Log
 All notable changes to this project will be documented in this file.
 
-## [0.5.0 - Gingold] - Mid August 2016
+## [0.5.1] - End of August
+
+### Added
+- Symmetric cell list
+- Verlet list
+- Full-Support for complex property on vector-dist (Serialization) + example
+
+### Fixed
+- Energy calculation in md example (double counting potential energy)
+
+### Changed
+
+## [0.5.0] - 15 August 2016
 
 ### Added
 - map communicate particles across processors mooving the information of all the particle map_list give the possibility to give a list of property to move from one to another processor

example/Vector/3_molecular_dynamic/main.cpp

@@ -218,11 +218,6 @@ double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd,
 		// Get the position of the particle p
 		Point<3,double> xp = vd.getPos(p);
 
-		// Reset the force
-		vd.template getProp<force>(p)[0] = 0.0;
-		vd.template getProp<force>(p)[1] = 0.0;
-		vd.template getProp<force>(p)[2] = 0.0;
-
 		// Get an iterator over the neighborhood of the particle p
 		auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
 
@@ -242,7 +237,7 @@ double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd,
 			double rn = norm2(xp - xq);
 
 			// potential energy (using pow is slower)
-			E += 4.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) );
+			E += 2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) );
 
 			// Next neighborhood
 			++Np;
@@ -286,8 +281,8 @@ int main(int argc, char* argv[])
 	//! \cond [constants] \endcond
 
 	double dt = 0.0005;
-	float r_cut = 0.3;
 	double sigma = 0.1;
+	double r_cut = 3.0*sigma;
 	double sigma12 = pow(sigma,12);
 	double sigma6 = pow(sigma,6);
 

example/Vector/3_molecular_dynamic/main_expr.cpp

@@ -15,7 +15,7 @@ struct ln_potential
 	{
 		double rn = norm2(xp - xq);
 
-		Point<2,double> E({4.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ),
+		Point<2,double> E({2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ),
 						 0.0});
 
 		return E;

openfpm_vcluster @ 0039a8c1

-Subproject commit a874640c9b69979cbf3117a4ecf9e471ae435b62
+Subproject commit 0039a8c1bc37029fdf99aaf633b1dfa64071fffc

src/Grid/grid_dist_id.hpp

@@ -519,15 +519,9 @@ class grid_dist_id
 		// check that the grid has valid size
 		check_size(g_sz);
 
-		// For a 5x5 grid you have 4x4 Cell (With the exception of periodic)
+		// get the size of the cell decomposer
 		size_t c_g[dim];
-		for (size_t i = 0 ; i < dim ; i++)
-		{
-			if (bc[i] == NON_PERIODIC)
-				c_g[i] = (g_sz[i]-1 > 0)?(g_sz[i]-1):1;
-			else
-				c_g[i] = g_sz[i];
-		}
+		getCellDecomposerPar<dim>(c_g,g_sz,bc);
 
 		// Initialize the cell decomposer
 		cd_sm.setDimensions(domain,c_g,0);

src/Grid/grid_dist_id_iterator_dec.hpp

@@ -145,6 +145,7 @@ class grid_dist_id_iterator_dec
 	 *
 	 * \param dec Decomposition
 	 * \param sz size of the grid
+	 * \param bc boundary conditions
 	 *
 	 */
 	grid_dist_id_iterator_dec(Decomposition & dec, const size_t (& sz)[Decomposition::dims])
@@ -152,7 +153,8 @@ class grid_dist_id_iterator_dec
 	{
 		// Initialize start and stop
 		start.zero();
-		for (size_t i = 0 ; i < Decomposition::dims ; i++) stop.set_d(i,sz[i]-1);
+		for (size_t i = 0 ; i < Decomposition::dims ; i++)
+			stop.set_d(i,sz[i]-1);
 
 		// From the decomposition construct gdb_ext
 		create_gdb_ext<Decomposition::dims,Decomposition>(gdb_ext,dec,sz,dec.getDomain(),spacing);
@@ -211,27 +213,6 @@ class grid_dist_id_iterator_dec
 		return *this;
 	}
 
-	/*! \brief Convert a g_dist_key_dx into a global key
-	 *
-	 * \see grid_dist_key_dx
-	 * \see grid_dist_iterator
-	 *
-	 * \return the global position in the grid
-	 *
-	 */
-/*	inline grid_key_dx<Decomposition::dims> getGKey(const grid_dist_key_dx<Decomposition::dims> & k)
-	{
-		// Get the sub-domain id
-		size_t sub_id = k.getSub();
-
-		grid_key_dx<Decomposition::dims> k_glob = k.getKey();
-
-		// shift
-		k_glob = k_glob + gdb_ext.get(sub_id).origin;
-
-		return k_glob;
-	}*/
-
 	/*! \brief Check if there is the next element
 	 *
 	 * \return true if there is the next, false otherwise

src/Grid/grid_dist_util.hpp

@@ -10,6 +10,25 @@
 
 #include "NN/CellList/CellDecomposer.hpp"
 
+/*! \brief get cellDecomposer parameters
+ *
+ * \tparam dim dimensionality
+ *
+ * \param c_g get the parameters of the cell decomposer
+ * \param g_sz global grid parameters
+ *
+ */
+template<unsigned int dim> void getCellDecomposerPar(size_t (& c_g)[dim], const size_t (& g_sz)[dim], const size_t (& bc)[dim])
+{
+	for (size_t i = 0 ; i < dim ; i++)
+	{
+		if (bc[i] == NON_PERIODIC)
+			c_g[i] = (g_sz[i]-1 > 0)?(g_sz[i]-1):1;
+		else
+			c_g[i] = g_sz[i];
+	}
+}
+
 /*! \brief
  *
  *
@@ -87,20 +106,21 @@ template<int dim, typename Decomposition> inline void create_gdb_ext(openfpm::ve
  * \param sz Global grid grid size
  * \param domain Domain where the grid is defined
  * \param spacing Define the spacing of the grid
+ * \param bc boundary conditions
  *
  */
 template<int dim, typename Decomposition> inline void create_gdb_ext(openfpm::vector<GBoxes<dim>> & gdb_ext, Decomposition & dec, const size_t (& sz)[dim], const Box<Decomposition::dims,typename Decomposition::stype> & domain, typename Decomposition::stype (& spacing)[dim])
 {
 	// Create the cell decomposer
-
 	CellDecomposer_sm<Decomposition::dims,typename Decomposition::stype, shift<Decomposition::dims,typename Decomposition::stype>> cd_sm;
 
-	size_t sz_cell[Decomposition::dims];
-	for (size_t i = 0 ; i < dim ; i++)
-		sz_cell[i] = sz[i] - 1;
+	size_t cdp[dim];
+
+	// Get the parameters to create a Cell-decomposer
+	getCellDecomposerPar<Decomposition::dims>(cdp,sz,dec.periodicity());
 
 	// Careful cd_sm require the number of cell
-	cd_sm.setDimensions(domain,sz_cell,0);
+	cd_sm.setDimensions(domain,cdp,0);
 
 	create_gdb_ext<dim,Decomposition>(gdb_ext,dec,cd_sm);
 

src/Vector/vector_dist.hpp

@@ -29,6 +29,7 @@
 #include "Vector/vector_dist_ofb.hpp"
 #include "Decomposition/CartDecomposition.hpp"
 #include "data_type/aggregate.hpp"
+#include "NN/VerletList/VerletList.hpp"
 
 #define V_SUB_UNIT_FACTOR 64
 
@@ -776,8 +777,7 @@ private:
 	 *
 	 * \return the processor bounding box
 	 */
-	inline Box<dim, St> cl_param_calculate(size_t (&div)[dim], St r_cut, const Ghost<dim, St> & enlarge)
-
+/*	inline Box<dim, St> cl_param_calculate(size_t (&div)[dim], St r_cut, const Ghost<dim, St> & enlarge)
 	{
 		// calculate the parameters of the cell list
 
@@ -795,7 +795,7 @@ private:
 			pbox.setHigh(i,pbox.getLow(i) + div[i]*r_cut);
 		}
 		return pbox;
-	}
+	}*/
 
 	/*! \brief Initialize the structures
 	 *
@@ -1456,8 +1456,11 @@ public:
 		// Division array
 		size_t div[dim];
 
+		// get the processor bounding box
+		Box<dim, St> pbox = dec.getProcessorBounds();
+
 		// Processor bounding box
-		auto pbox = cl_param_calculate(div, r_cut, enlarge);
+		cl_param_calculate(pbox, div, r_cut, enlarge);
 
 		cell_list.Initialize(pbox, div);
 
@@ -1488,8 +1491,11 @@ public:
 		// Division array
 		size_t div[dim];
 
+		// get the processor bounding box
+		Box<dim, St> pbox = dec.getProcessorBounds();
+
 		// Processor bounding box
-		auto pbox = cl_param_calculate(div, r_cut, enlarge);
+		cl_param_calculate(pbox,div, r_cut, enlarge);
 
 		cell_list.Initialize(pbox, div, g_m);
 
@@ -1498,13 +1504,39 @@ public:
 		return cell_list;
 	}
 
+
+	/*! \brief for each particle get the verlet list
+	 *
+	 * \param r_cut cut-off radius
+	 *
+	 */
+	VerletList<dim,St,FAST,shift<dim,St> > getVerlet(St r_cut)
+	{
+		VerletList<dim,St,FAST,shift<dim,St>> ver;
+
+		// Division array
+		size_t div[dim];
+
+		// get the processor bounding box
+		Box<dim, St> bt = dec.getProcessorBounds();
+
+		// Calculate the divisions for the Cell-lists
+		cl_param_calculate(bt,div,r_cut,Ghost<dim,St>(0.0));
+
+		ver.Initialize(bt,r_cut,v_pos,g_m);
+
+		return ver;
+	}
+
 	/*! \brief for each particle get the verlet list
 	 *
 	 * \param verlet output verlet list for each particle
 	 * \param r_cut cut-off radius
 	 *
+	 * \deprecated
+	 *
 	 */
-	void getVerlet(openfpm::vector<openfpm::vector<size_t>> & verlet, St r_cut)
+	void getVerletDeprecated(openfpm::vector<openfpm::vector<size_t>> & verlet, St r_cut)
 	{
 		// resize verlet to store the number of particles
 		verlet.resize(size_local());
@@ -1722,25 +1754,15 @@ public:
 	 */
 	inline grid_dist_id_iterator_dec<Decomposition> getGridIterator(const size_t (&sz)[dim])
 	{
-		size_t sz_g[dim];
 		grid_key_dx<dim> start;
 		grid_key_dx<dim> stop;
 		for (size_t i = 0; i < dim; i++)
 		{
 			start.set_d(i, 0);
-			if (dec.periodicity(i) == PERIODIC)
-			{
-				sz_g[i] = sz[i];
-				stop.set_d(i, sz_g[i] - 2);
-			}
-			else
-			{
-				sz_g[i] = sz[i];
-				stop.set_d(i, sz_g[i] - 1);
-			}
+			stop.set_d(i, sz[i] - 1);
 		}
 
-		grid_dist_id_iterator_dec<Decomposition> it_dec(dec, sz_g, start, stop);
+		grid_dist_id_iterator_dec<Decomposition> it_dec(dec, sz, start, stop);
 		return it_dec;
 	}
 

src/Vector/vector_dist_unit_test.hpp

@@ -1086,6 +1086,71 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_interacting_particles )
 	}
 }
 
+BOOST_AUTO_TEST_CASE( vector_dist_grid_iterator )
+{
+	long int k = 64*64*64*create_vcluster().getProcessingUnits();
+	k = std::pow(k, 1/3.);
+
+	long int big_step = k / 30;
+	big_step = (big_step == 0)?1:big_step;
+	long int small_step = 21;
+
+	print_test( "Testing vector grid iterator list k<=",k);
+
+	// 3D test
+	for ( ; k > 8*big_step ; k-= (k > 2*big_step)?big_step:small_step )
+	{
+		Vcluster & v_cl = create_vcluster();
+
+		const size_t Ng = k;
+
+		// we create a 128x128x128 Grid iterator
+		size_t sz[3] = {Ng,Ng,Ng};
+
+		Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+		// Boundary conditions
+		size_t bc[3]={NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
+
+		// ghost
+		Ghost<3,float> ghost(1.0/(Ng-2));
+
+		// Distributed vector
+		vector_dist<3,float, Point_test<float>, CartDecomposition<3,float> > vd(0,box,bc,ghost);
+
+		// Put particles on a grid creating a Grid iterator
+		auto it = vd.getGridIterator(sz);
+
+		while (it.isNext())
+		{
+			vd.add();
+
+			auto key = it.get();
+
+			vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
+			vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
+			vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
+
+			++it;
+		}
+
+		BOOST_REQUIRE_EQUAL(it.getSpacing(0),1.0f/(Ng-1));
+		BOOST_REQUIRE_EQUAL(it.getSpacing(1),1.0f/(Ng-1));
+		BOOST_REQUIRE_EQUAL(it.getSpacing(2),1.0f/(Ng-1));
+
+		// distribute particles and sync ghost
+		vd.map();
+
+
+		// Check that the sum of all the particles is the grid size
+		size_t total = vd.size_local();
+		v_cl.sum(total);
+		v_cl.execute();
+
+		BOOST_REQUIRE_EQUAL(total,(Ng) * (Ng) * (Ng));
+	}
+}
+
 BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
 {
 	long int k = 64*64*64*create_vcluster().getProcessingUnits();
@@ -1134,6 +1199,10 @@ BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
 			++it;
 		}
 
+		BOOST_REQUIRE_EQUAL(it.getSpacing(0),1.0f/Ng);
+		BOOST_REQUIRE_EQUAL(it.getSpacing(1),1.0f/Ng);
+		BOOST_REQUIRE_EQUAL(it.getSpacing(2),1.0f/Ng);
+
 		// distribute particles and sync ghost
 		vd.map();
 
@@ -1142,7 +1211,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
 		v_cl.sum(total);
 		v_cl.execute();
 
-		BOOST_REQUIRE_EQUAL(total,(Ng-1) * (Ng-1) * (Ng-1));
+		BOOST_REQUIRE_EQUAL(total,(Ng) * (Ng) * (Ng));
 
 		vd.ghost_get<0>();
 
@@ -1162,9 +1231,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
 
 		// Create a verlet list for each particle
 
-		openfpm::vector<openfpm::vector<size_t>> verlet;
-
-		vd.getVerlet(verlet,third_dist);
+		VerletList<3,float,FAST,shift<3,float>> verlet = vd.getVerlet(third_dist);
 
 		bool correct = true;
 
@@ -1179,11 +1246,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
 			Point<3,float> p = vd.getPos(i);
 
 			// for each neighborhood particle
-			for (size_t j = 0 ; j < verlet.get(i).size() ; j++)
+			for (size_t j = 0 ; j < verlet.getNNPart(i) ; j++)
 			{
-				auto & NN = verlet.get(i);
-
-				Point<3,float> q = vd.getPos(NN.get(j));
+				Point<3,float> q = vd.getPos(verlet.get(i,j));
 
 				float dist = p.distance(q);
 
@@ -1195,7 +1260,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
 					third_NN++;
 			}
 
-			correct &= (first_NN == 6);
+			correct &= (first_NN == 7);
 			correct &= (second_NN == 12);
 			correct &= (third_NN == 8);
 		}
@@ -1324,323 +1389,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_map_list )
 	}
 }
 
-///////////////////////// test hilb ///////////////////////////////
-
-BOOST_AUTO_TEST_CASE( vector_dist_reorder_2d_test )
-{
-	Vcluster & v_cl = create_vcluster();
-
-    // set the seed
-	// create the random generator engine
-	std::srand(v_cl.getProcessUnitID());
-    std::default_random_engine eg;
-    std::uniform_real_distribution<float> ud(0.0f, 1.0f);
-
-    long int k = 524288 * v_cl.getProcessingUnits();
-
-	long int big_step = k / 4;
-	big_step = (big_step == 0)?1:big_step;
-
-	print_test_v( "Testing 2D vector with hilbert curve reordering k<=",k);
-
-	// 2D test
-	for ( ; k >= 2 ; k-= decrement(k,big_step) )
-	{
-		BOOST_TEST_CHECKPOINT( "Testing 2D vector with hilbert curve reordering k=" << k );
-
-		Box<2,float> box({0.0,0.0},{1.0,1.0});
-
-		// Boundary conditions
-		size_t bc[2]={NON_PERIODIC,NON_PERIODIC};
-
-		vector_dist<2,float, Point_test<float>, CartDecomposition<2,float> > vd(k,box,bc,Ghost<2,float>(0.0));
-
-		auto it = vd.getIterator();
-
-		while (it.isNext())
-		{
-			auto key = it.get();
-
-			vd.getPos(key)[0] = ud(eg);
-			vd.getPos(key)[1] = ud(eg);
-
-			++it;
-		}
-
-		vd.map();
-
-		// Create first cell list
-
-		auto NN1 = vd.getCellList(0.01);
-
-		//An order of a curve
-		int32_t m = 6;
-
-		//Reorder a vector
-		vd.reorder(m);
-
-		// Create second cell list
-		auto NN2 = vd.getCellList(0.01);
-
-		//Check equality of cell sizes
-		for (size_t i = 0 ; i < NN1.getGrid().size() ; i++)
-		{
-			size_t n1 = NN1.getNelements(i);
-			size_t n2 = NN2.getNelements(i);
-
-			BOOST_REQUIRE_EQUAL(n1,n2);
-		}
-	}
-}
-
-BOOST_AUTO_TEST_CASE( vector_dist_cl_random_vs_hilb_forces_test )
-{
-	Vcluster & v_cl = create_vcluster();
-
-	if (v_cl.getProcessingUnits() > 32)
-		return;
-
-	///////////////////// INPUT DATA //////////////////////
-
-	// Dimensionality of the space
-	const size_t dim = 3;
-	// Cut-off radiuses. Can be put different number of values
-	openfpm::vector<float> cl_r_cutoff {0.05};
-	// The starting amount of particles (remember that this number is multiplied by number of processors you use for testing)
-	size_t cl_k_start = 10000;
-	// The lower threshold for number of particles
-	size_t cl_k_min = 1000;
-	// Ghost part of distributed vector
-	double ghost_part = 0.05;
-
-	///////////////////////////////////////////////////////
-
-	//For different r_cut
-	for (size_t r = 0; r < cl_r_cutoff.size(); r++ )
-	{
-		//Cut-off radius
-		float r_cut = cl_r_cutoff.get(r);
-
-		//Number of particles
-		size_t k = cl_k_start * v_cl.getProcessingUnits();
-
-		std::string str("Testing " + std::to_string(dim) + "D vector's forces (random vs hilb celllist) k<=");
-
-		vector_dist_test::print_test_v(str,k);
-
-		//For different number of particles
-		for (size_t k_int = k ; k_int >= cl_k_min ; k_int/=2 )
-		{
-			BOOST_TEST_CHECKPOINT( "Testing " << dim << "D vector's forces (random vs hilb celllist) k<=" << k_int );
-
-			Box<dim,float> box;
-
-			for (size_t i = 0; i < dim; i++)
-			{
-				box.setLow(i,0.0);
-				box.setHigh(i,1.0);
-			}
-
-			// Boundary conditions
-			size_t bc[dim];
-
-			for (size_t i = 0; i < dim; i++)
-				bc[i] = PERIODIC;
-
-			vector_dist<dim,float, aggregate<float[dim]>, CartDecomposition<dim,float> > vd(k_int,box,bc,Ghost<dim,float>(ghost_part));
-
-			vector_dist<dim,float, aggregate<float[dim]>, CartDecomposition<dim,float> > vd2(k_int,box,bc,Ghost<dim,float>(ghost_part));
-
-			// Initialize dist vectors
-			vd_initialize_double<dim>(vd, vd2, v_cl, k_int);
-
-			vd.template ghost_get<0>();
-			vd2.template ghost_get<0>();
-
-			//Get a cell list
-
-			auto NN = vd.getCellList(r_cut);
-
-			//Calculate forces
-
-			calc_forces<dim>(NN,vd,r_cut);
-
-			//Get a cell list hilb
-
-			auto NN_hilb = vd2.getCellList_hilb(r_cut);
-
-			//Calculate forces
-			calc_forces_hilb<dim>(NN_hilb,vd2,r_cut);
-
-			// Calculate average
-			size_t count = 1;
-			Point<dim,float> avg;
-			for (size_t i = 0 ; i < dim ; i++)	{avg.get(i) = 0.0;}
-
-			auto it_v2 = vd.getIterator();
-			while (it_v2.isNext())
-			{
-				//key
-				vect_dist_key_dx key = it_v2.get();
-
-				for (size_t i = 0; i < dim; i++)
-					avg.get(i) += fabs(vd.template getProp<0>(key)[i]);
-
-				++count;
-				++it_v2;
-			}
-
-			for (size_t i = 0 ; i < dim ; i++)	{avg.get(i) /= count;}
-
-			auto it_v = vd.getIterator();
-			while (it_v.isNext())
-			{
-				//key
-				vect_dist_key_dx key = it_v.get();
-
-				for (size_t i = 0; i < dim; i++)