Workinh CRS with debug garbage

src/Decomposition/CartDecomposition.hpp

@@ -36,9 +36,13 @@
 #include "util/se_util.hpp"
 #include "util/mathutil.hpp"
 #include "CartDecomposition_ext.hpp"
+#include "data_type/aggregate.hpp"
+#include "Domain_NN_calculator_cart.hpp"
 
 #define CARTDEC_ERROR 2000lu
 
+#define COMPUTE_SKIN_SUB 1
+
 /**
  * \brief This class decompose a space into sub-sub-domains and distribute them across processors
  *
@@ -81,7 +85,7 @@
  */
 
 template<unsigned int dim, typename T, typename Memory, typename Distribution>
-class CartDecomposition: public ie_loc_ghost<dim, T>, public nn_prcs<dim, T>, public ie_ghost<dim, T>
+class CartDecomposition: public ie_loc_ghost<dim, T>, public nn_prcs<dim, T>, public ie_ghost<dim, T>, public domain_nn_calculator_cart<dim>
 {
 
 public:
@@ -153,6 +157,12 @@ protected:
 	//! Boundary condition info
 	size_t bc[dim];
 
+	//! Processor domain bounding box
+	::Box<dim,size_t> proc_box;
+
+	//! set of Boxes produced by the decomposition optimizer
+	openfpm::vector<::Box<dim, size_t>> loc_box;
+
 	/*! \brief It convert the box from the domain decomposition into sub-domain
 	 *
 	 * The decomposition box from the domain-decomposition contain the box in integer
@@ -205,9 +215,10 @@ public:
 	 *
 	 * \param v_cl Virtual cluster, used internally for communications
 	 * \param bc boundary conditions
+	 * \param opt option (one option is to construct)
 	 *
 	 */
-	void createSubdomains(Vcluster & v_cl, const size_t (& bc)[dim])
+	void createSubdomains(Vcluster & v_cl, const size_t (& bc)[dim], size_t opt = 0)
 	{
 #ifdef SE_CLASS1
 		if (&v_cl == NULL)
@@ -237,9 +248,6 @@ public:
 		// And reducing Ghost over-stress
 		dec_optimizer<dim, Graph_CSR<nm_v, nm_e>> d_o(dist.getGraph(), gr.getSize());
 
-		// set of Boxes produced by the decomposition optimizer
-		openfpm::vector<::Box<dim, size_t>> loc_box;
-
 		// Ghost
 		Ghost<dim,long int> ghe;
 
@@ -259,7 +267,10 @@ public:
 
 		// Initialize ss_box and bbox
 		if (loc_box.size() >= 0)
+		{
 			bbox = convertDecBoxIntoSubDomain(loc_box.get(0));
+			proc_box = loc_box.get(0);
+		}
 
 		// convert into sub-domain
 		for (size_t s = 0; s < loc_box.size(); s++)
@@ -271,6 +282,7 @@ public:
 
 			// Calculate the bound box
 			bbox.enclose(sub_d);
+			proc_box.enclose(loc_box.get(s));
 
 			// Create the smallest box contained in all sub-domain
 			ss_box.contained(sub_d);
@@ -362,7 +374,7 @@ public:
 		}
 	}
 
-	/*! \brief Create the subspaces that decompose your domain
+	/*! \brief Create the sub-domain that decompose your domain
 	 *
 	 */
 	void CreateSubspaces()
@@ -392,7 +404,7 @@ public:
 			//! add the space box
 			sub_domains.add(tmp);
 
-			// add the iterator
+			// Next sub-domain
 			++gk_it;
 		}
 	}
@@ -531,8 +543,8 @@ public:
 	 * \param v_cl Virtual cluster, used internally to handle or pipeline communication
 	 *
 	 */
-	CartDecomposition(Vcluster & v_cl) :
-			nn_prcs<dim, T>(v_cl), v_cl(v_cl), dist(v_cl),ref_cnt(0)
+	CartDecomposition(Vcluster & v_cl)
+	:nn_prcs<dim, T>(v_cl), v_cl(v_cl), dist(v_cl),ref_cnt(0)
 	{
 		// Reset the box to zero
 		bbox.zero();
@@ -1065,10 +1077,6 @@ public:
 			{
 				std::cout << std::setprecision(3) << unbalance << "\n";
 			}
-
-//			write(v_cl.getProcessUnitID() + "_"+ std::to_string(n_step) + "_AAAAAA");
-
-//			n_step++;
 		}
 
 		if (dlb.rebalanceNeeded())
@@ -1268,6 +1276,28 @@ public:
 		return processorID<Mem>(pt) == v_cl.getProcessUnitID();
 	}
 
+	/*! \brief Get the domain Cells
+	 *
+	 * \param shift Shifting point
+	 * \param gs grid extension
+	 *
+	 */
+	openfpm::vector<size_t> & getDomainCells(grid_key_dx<dim> & shift, grid_key_dx<dim> & cell_shift, grid_sm<dim,void> & gs)
+	{
+		return domain_nn_calculator_cart<dim>::getDomainCells(shift,cell_shift,gs,proc_box,loc_box);
+	}
+
+	/*! \brief Get the domain anomalous cells
+	 *
+	 * \param shift Shifting point
+	 * \param gs grid extension
+	 *
+	 */
+	openfpm::vector<subsub_lin<dim>> & getAnomDomainCells(grid_key_dx<dim> & shift, grid_key_dx<dim> & cell_shift, grid_sm<dim,void> & gs)
+	{
+		return domain_nn_calculator_cart<dim>::getAnomDomainCells(shift,cell_shift,gs,proc_box,loc_box);
+	}
+
 	/*! \brief Check if the particle is local considering boundary conditions
 	 *
 	 * \warning if the particle id outside the domain and non periodic boundary the result

src/Graph/dist_map_graph_unit_test.hpp

@@ -58,9 +58,9 @@ struct vx
 
 	template<unsigned int dim, typename Mem> inline vx & operator=(const encapc<dim, vx, Mem> & p)
 	{
-		boost::fusion::at_c<0>(data)[0] = p.get<0>()[0];
-		boost::fusion::at_c<0>(data)[1] = p.get<0>()[1];
-		boost::fusion::at_c<0>(data)[2] = p.get<0>()[2];
+		boost::fusion::at_c<0>(data)[0] = p.template get<0>()[0];
+		boost::fusion::at_c<0>(data)[1] = p.template get<0>()[1];
+		boost::fusion::at_c<0>(data)[2] = p.template get<0>()[2];
 
 		return *this;
 	}

src/Vector/vector_dist.hpp

@@ -113,6 +113,19 @@ private:
 		g_m = p_np;
 	}
 
+	/*! \brief Checl if the parameters describe a valid vector. In case it does not report an error
+	 *
+	 * \param box Box to check
+	 *
+	 */
+	void check_parameters(Box<dim,St> & box)
+	{
+		// if the box is not valid return an error
+		if (box.isValid() == false)
+			std::cerr << __FILE__ << ":" << __LINE__ << "  Error the domain is not valid " << box.toString() << std::endl;
+
+	}
+
 public:
 
 	//! space type
@@ -223,6 +236,8 @@ public:
 		check_new(this,8,VECTOR_DIST_EVENT,4);
 #endif
 
+		check_parameters(box);
+
 		init_structures(np);
 		this->init_decomposition(box,bc,g,opt);
 	}
@@ -1157,6 +1172,60 @@ public:
 
 		return sz;
 	}
+
+	/*! \brief Get a special particle iterator able to iterate across particles using
+	 *         symmetric crossing scheme
+	 *
+	 * \param NN Cell-List neighborhood
+	 *
+	 */
+	template<typename cli> ParticleItCRS_Cells<dim,cli> getParticleIteratorCRS(cli & NN)
+	{
+		// Shift
+		grid_key_dx<dim> cell_shift = NN.getShift();
+
+		// Shift
+		grid_key_dx<dim> shift = NN.getShift();
+
+		// Add padding
+		for (size_t i = 0 ; i < dim ; i++)
+			shift.set_d(i,shift.get(i) + NN.getPadding(i));
+
+		grid_sm<dim,void> gs = NN.getInternalGrid();
+
+		// First we check that
+		return ParticleItCRS_Cells<dim,cli>(NN,getDecomposition().getDomainCells(shift,cell_shift,gs),getDecomposition().getAnomDomainCells(shift,cell_shift,gs),NN.getNNc_sym());
+	}
+
+	/*! \brief Return from which cell we have to start in case of CRS interation
+	 *         scheme
+	 *
+	 * \param NN cell-list
+	 *
+	 * \return The starting cell point
+	 *
+	 */
+	template<typename Celllist> grid_key_dx<dim> getCRSStart(Celllist & NN)
+	{
+		return NN.getStartDomainCell();
+	}
+
+	/*! \brief Return from which cell we have to stop in case of CRS interation
+	 *         scheme
+	 *
+	 * \param NN cell-list
+	 *
+	 * \return The stop cell point
+	 *
+	 */
+	template<typename Celllist> grid_key_dx<dim> getCRSStop(Celllist & NN)
+	{
+		grid_key_dx<dim> key = NN.getStopDomainCell();
+
+		for (size_t i = 0 ; i < dim ; i++)
+			key.set_d(i,key.get(i) + 1);
+		return key;
+	}
 };
 
 

src/Vector/vector_dist_cell_list_tests.hpp

@@ -533,6 +533,261 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
 	BOOST_REQUIRE_EQUAL(ret,true);
 }
 
+BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
+{
+	Vcluster & v_cl = create_vcluster();
+
+	if (v_cl.getProcessingUnits() > 24)
+		return;
+
+	float L = 1000.0;
+
+    // set the seed
+	// create the random generator engine
+    std::default_random_engine eg(1132312*v_cl.getProcessUnitID());
+    std::uniform_real_distribution<float> ud(-L,L);
+
+    long int k = 4096 * v_cl.getProcessingUnits();
+
+	long int big_step = k / 4;
+	big_step = (big_step == 0)?1:big_step;
+
+	print_test("Testing 3D periodic vector symmetric cell-list k=",k);
+	BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list k=" << k );
+
+	Box<3,float> box({-L,-L,-L},{L,L,L});
+
+	// Boundary conditions
+	size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
+
+	float r_cut = 100.0;
+
+	// ghost
+	Ghost<3,float> ghost(r_cut);
+
+	// Point and global id
+	struct point_and_gid
+	{
+		size_t id;
+		Point<3,float> xq;
+
+		bool operator<(const struct point_and_gid & pag) const
+		{
+			return (id < pag.id);
+		}
+	};
+
+	typedef  aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
+
+	// Distributed vector
+	vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_GHOST);
+	size_t start = vd.init_size_accum(k);
+
+	auto it = vd.getIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		vd.getPos(key)[0] = ud(eg);
+		vd.getPos(key)[1] = ud(eg);
+		vd.getPos(key)[2] = ud(eg);
+
+		// Fill some properties randomly
+
+		vd.getProp<0>(key) = 0;
+		vd.getProp<1>(key) = 0;
+		vd.getProp<2>(key) = key.getKey() + start;
+
+		++it;
+	}
+
+	vd.map();
+
+	// sync the ghost
+	vd.ghost_get<0,2>();
+
+	auto NN = vd.getCellList(r_cut);
+	auto p_it = vd.getDomainIterator();
+
+	while (p_it.isNext())
+	{
+		auto p = p_it.get();
+
+		Point<3,float> xp = vd.getPos(p);
+
+		auto Np = NN.getNNIterator(NN.getCell(vd.getPos(p)));
+
+		while (Np.isNext())
+		{
+			auto q = Np.get();
+
+			if (p.getKey() == q)
+			{
+				++Np;
+				continue;
+			}
+
+			// repulsive
+
+			Point<3,float> xq = vd.getPos(q);
+			Point<3,float> f = (xp - xq);
+
+			float distance = f.norm();
+
+			// Particle should be inside 2 * r_cut range
+
+			if (distance < r_cut )
+			{
+				vd.getProp<0>(p)++;
+				vd.getProp<3>(p).add();
+				vd.getProp<3>(p).last().xq = xq;
+				vd.getProp<3>(p).last().id = vd.getProp<2>(q);
+			}
+
+			++Np;
+		}
+
+		++p_it;
+	}
+
+	// We now try symmetric  Cell-list
+
+	auto NN2 = vd.getCellListSym(r_cut);
+
+	int debug = 0;
+	debug++;
+
+	auto debug_it = vd.getDomainAndGhostIterator();
+
+	while (debug_it.isNext())
+	{
+		auto key = debug_it.get();
+
+		if (vd.getProp<2>(key) == 1698)
+		{
+			int debug = 0;
+			debug++;
+		}
+
+		++debug_it;
+	}
+
+	vd.write("debug_decomp_crs_part");
+	vd.getDecomposition().write("debug_decomp_crs");
+	std::cout << "NN2:   " << NN2.getCell(vd.getPos(4853)) << std::endl;
+	std::cout << "NN2 Grid:   " << NN2.getCellGrid(vd.getPos(4853)).to_string() << std::endl;
+
+	int debug_cnt = 0;
+
+	// In case of CRS we have to iterate particles within some cells
+	// here we define whichone
+	auto p_it2 = vd.getParticleIteratorCRS(NN2);
+
+	// For each particle
+	while (p_it2.isNext())
+	{
+		auto p = p_it2.get();
+
+		debug_cnt++;
+		Point<3,float> xp = vd.getPos(p);
+
+		if (debug_cnt == 3762  && v_cl.getProcessUnitID() == 0)
+		{
+			std::cout << " P local " << p << "      " << vd.getProp<2>(p) << std::endl;
+			int debug = 0;
+			debug++;
+		}
+
+		auto Np = p_it2.getNNIteratorCSR(vd.getPosVector());
+
+		while (Np.isNext())
+		{
+			auto q = Np.get();
+
+			if (p == q)
+			{
+				++Np;
+				continue;
+			}
+
+			// repulsive
+
+			Point<3,float> xq = vd.getPos(q);
+			Point<3,float> f = (xp - xq);
+
+			float distance = f.norm();
+
+			// Particle should be inside r_cut range
+
+			if (distance < r_cut )
+			{
+				vd.getProp<1>(p)++;
+				vd.getProp<1>(q)++;
+
+				vd.getProp<4>(p).add();
+				vd.getProp<4>(q).add();
+
+				vd.getProp<4>(p).last().xq = xq;
+				vd.getProp<4>(q).last().xq = xp;
+				vd.getProp<4>(p).last().id = vd.getProp<2>(q);
+				vd.getProp<4>(q).last().id = vd.getProp<2>(p);
+
+				if (v_cl.getProcessUnitID() == 0 && vd.getProp<2>(p) == 1698)
+				{
+					std::cerr << "NN " << vd.getProp<2>(q) << std::endl;
+				}
+
+				if (v_cl.getProcessUnitID() == 0 && vd.getProp<2>(q) == 1689)
+				{
+					std::cerr << "NN " << vd.getProp<2>(p) << std::endl;
+				}
+			}
+
+			++Np;
+		}
+
+		++p_it2;
+	}
+
+	vd.ghost_put<add_,1>();
+	vd.ghost_put<merge_,4>();
+
+	auto p_it3 = vd.getDomainIterator();
+
+	bool ret = true;
+	while (p_it3.isNext())
+	{
+		auto p = p_it3.get();
+
+		ret &= vd.getProp<1>(p) == vd.getProp<0>(p);
+
+		vd.getProp<3>(p).sort();
+		vd.getProp<4>(p).sort();
+
+		ret &= vd.getProp<3>(p).size() == vd.getProp<4>(p).size();
+
+		if (v_cl.getProcessUnitID() == 0)
+		{
+			std::cerr << "Particle: " << p.getKey()  <<  "   Position: " << Point<3,float>(vd.getPos(p)).toString() << std::endl;
+
+			for (size_t i = 0 ; i < vd.getProp<4>(p).size() ; i++)
+			{
+				std::cerr << "POSITION: " << vd.getProp<3>(p).get(i).xq.toString() << std::endl;
+
+				std::cerr << "ID nn " << vd.getProp<3>(p).get(i).id << "     " << vd.getProp<4>(p).get(i).id << std::endl;
+				ret &= vd.getProp<3>(p).get(i).id == vd.getProp<4>(p).get(i).id;
+			}
+		}
+
+		if (ret == false)
+			break;
+
+		++p_it3;
+	}
+
+	BOOST_REQUIRE_EQUAL(ret,true);
+}
 
 BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list )
 {
@@ -956,4 +1211,5 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list_no_bottom )
 	}
 }
 
+
 #endif /* SRC_VECTOR_VECTOR_DIST_CELL_LIST_TESTS_HPP_ */

src/Vector/vector_dist_comm.hpp

@@ -919,9 +919,6 @@ public:
 		// Processor communication size
 		openfpm::vector<size_t> prc_sz(v_cl.getProcessingUnits());
 
-		// It contain the list of the processors this processor should to communicate with
-		openfpm::vector<size_t> p_list;
-
 		// map completely reset the ghost part
 		v_pos.resize(g_m);
 		v_prp.resize(g_m);
@@ -981,9 +978,6 @@ public:
 		// Processor communication size
 		openfpm::vector<size_t> prc_sz(v_cl.getProcessingUnits());
 
-		// It contain the list of the processors this processor should to communicate with
-		openfpm::vector<size_t> p_list;
-
 		// map completely reset the ghost part
 		v_pos.resize(g_m);
 		v_prp.resize(g_m);