diff --git a/openfpm_numerics b/openfpm_numerics index 71fc012b7e74ece8ac880dddbc9b9aabdb6d3568..55add63fdd24344834dc27de3028ed0c09c0e6d2 160000 --- a/openfpm_numerics +++ b/openfpm_numerics @@ -1 +1 @@ -Subproject commit 71fc012b7e74ece8ac880dddbc9b9aabdb6d3568 +Subproject commit 55add63fdd24344834dc27de3028ed0c09c0e6d2 diff --git a/src/Grid/grid_dist_id.hpp b/src/Grid/grid_dist_id.hpp index 7a0f96cc0590f8ff38fd9b77c83530d56e28d34f..ca7f207f7cb052498246c066b885cd582c8275d8 100644 --- a/src/Grid/grid_dist_id.hpp +++ b/src/Grid/grid_dist_id.hpp @@ -1109,7 +1109,7 @@ public: struct i_box_id { //! Box - ::Box<dim,size_t> box; + ::Box<dim,long int> box; //! id size_t g_id; @@ -1128,7 +1128,7 @@ public: struct i_lbox_id { //! Box - ::Box<dim,size_t> box; + ::Box<dim,long int> box; //! sub-domain id size_t sub; @@ -1144,10 +1144,10 @@ public: struct e_box_id { //! Box defining the external ghost box in global coordinates - ::Box<dim,size_t> g_e_box; + ::Box<dim,long int> g_e_box; //! Box defining the external ghost box in local coordinates - ::Box<dim,size_t> l_e_box; + ::Box<dim,long int> l_e_box; //! Sector position of the external ghost comb<dim> cmb; @@ -1166,7 +1166,7 @@ public: struct e_lbox_id { //! Box defining the external ghost box in local coordinates - ::Box<dim,size_t> box; + ::Box<dim,long int> box; //! Sector position of the local external ghost box comb<dim> cmb; @@ -1282,7 +1282,7 @@ public: // And linked sub-domain size_t sub_id = ig_box.get(i).bid.get(j).sub; // Internal ghost box - Box<dim,size_t> g_ig_box = ig_box.get(i).bid.get(j).box; + Box<dim,long int> g_ig_box = ig_box.get(i).bid.get(j).box; if (g_ig_box.isValid() == false) continue; @@ -1521,6 +1521,37 @@ public: #endif } + /*! \brief It print the internal ghost boxes and external ghost boxes in global unit + * + * + */ + void debugPrint() + { + std::cout << "-------- External Ghost boxes ---------- " << std::endl; + + for (size_t i = 0 ; i < eg_box.size() ; i++) + { + std::cout << "Processor: " << eg_box.get(i).prc << " Boxes:" << std::endl; + + for (size_t j = 0; j < eg_box.get(i).bid.size() ; j++) + { + std::cout << " Box: " << eg_box.get(i).bid.get(j).g_e_box.toString() << " Id: " << eg_box.get(i).bid.get(j).g_id << std::endl; + } + } + + std::cout << "-------- Internal Ghost boxes ---------- " << std::endl; + + for (size_t i = 0 ; i < ig_box.size() ; i++) + { + std::cout << "Processor: " << ig_box.get(i).prc << " Boxes:" << std::endl; + + for (size_t j = 0 ; j < ig_box.get(i).bid.size() ; j++) + { + std::cout << " Box: " << ig_box.get(i).bid.get(j).box.toString() << " Id: " << ig_box.get(i).bid.get(j).g_id << std::endl; + } + } + } + // Define friend classes friend grid_dist_id<dim,St,T,typename Decomposition::extended_type,Memory,device_grid>; diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp index 0e28d613fe032477fd32e14b822b3a4478afd785..8e79e5be8cbdb7aa4098f19247de670fc0e2136c 100644 --- a/src/Vector/vector_dist.hpp +++ b/src/Vector/vector_dist.hpp @@ -137,6 +137,15 @@ private: openfpm::vector<prop, PreAllocHeapMemory<2>, typename memory_traits_lin<prop>::type, memory_traits_lin, openfpm::grow_policy_identity> prp; }; + template <typename sel_prop> + struct pos_prop_sel + { + //! position vector + openfpm::vector<Point<dim, St>, PreAllocHeapMemory<2>, typename memory_traits_lin<Point<dim, St>>::type, memory_traits_lin, openfpm::grow_policy_identity> pos; + //! properties vector + openfpm::vector<sel_prop, PreAllocHeapMemory<2>, typename memory_traits_lin<sel_prop>::type, memory_traits_lin, openfpm::grow_policy_identity> prp; + }; + /*! \brief Label particles for mappings * * \param lbl_p Particle labeled @@ -500,6 +509,63 @@ private: } } + /*! \brief allocate and fill the send buffer for the map function + * + * \param prc_r List of processor rank involved in the send + * \param prc_r_sz For each processor in the list the size of the message to send + * \param pb send buffer + * + */ + template<typename prp_object,int ... prp> void fill_send_map_buf_list(openfpm::vector<size_t> & prc_r, openfpm::vector<size_t> & prc_sz_r, openfpm::vector<pos_prop_sel<prp_object>> & pb) + { + pb.resize(prc_r.size()); + + for (size_t i = 0; i < prc_r.size(); i++) + { + // Create the size required to store the particles position and properties to communicate + size_t s1 = openfpm::vector<Point<dim, St>, HeapMemory, typename memory_traits_lin<Point<dim, St>>::type, memory_traits_lin, openfpm::grow_policy_identity>::calculateMem(prc_sz_r.get(i), 0); + size_t s2 = openfpm::vector<prp_object, HeapMemory, typename memory_traits_lin<prp_object>::type, memory_traits_lin, openfpm::grow_policy_identity>::calculateMem(prc_sz_r.get(i), 0); + + // Preallocate the memory + size_t sz[2] = { s1, s2 }; + PreAllocHeapMemory<2> * mem = new PreAllocHeapMemory<2>(sz); + + // Set the memory allocator + pb.get(i).pos.setMemory(*mem); + pb.get(i).prp.setMemory(*mem); + + // set the size and allocate, using mem warant that pos and prp is contiguous + pb.get(i).pos.resize(prc_sz_r.get(i)); + pb.get(i).prp.resize(prc_sz_r.get(i)); + } + + // Run through all the particles and fill the sending buffer + + for (size_t i = 0; i < opart.size(); i++) + { + auto it = opart.get(i).getIterator(); + size_t lbl = p_map_req.get(i); + + while (it.isNext()) + { + size_t key = it.get(); + size_t id = opart.get(i).get(key); + + pb.get(lbl).pos.set(key, v_pos.get(id)); + + // source object type + typedef encapc<1, prop, typename openfpm::vector<prop>::layout_type> encap_src; + // destination object type + typedef encapc<1, prp_object, typename openfpm::vector<prp_object>::layout_type> encap_dst; + + // Copy only the selected properties + object_si_d<encap_src, encap_dst, OBJ_ENCAP, prp...>(v_prp.get(id), pb.get(lbl).prp.get(key)); + + ++it; + } + } + } + /*! \brief This function process the receiced data for the properties and populate the ghost * * \tparam send_vector type used to send data @@ -620,6 +686,68 @@ private: v_prp.remove(out_part, o_p_id); } + /*! \brief Process the received particles + * + * \param list of the out-going particles + * + */ + template<typename prp_object , int ... prp> void process_received_map_list(openfpm::vector<size_t> & out_part) + { + size_t o_p_id = 0; + + for (size_t i = 0; i < recv_mem_gm.size(); i++) + { + // Get the number of elements + + size_t n_ele = recv_mem_gm.get(i).size() / (sizeof(Point<dim, St> ) + sizeof(prop)); + + // Pointer of the received positions for each near processor + void * ptr_pos = (unsigned char *) recv_mem_gm.get(i).getPointer(); + // Pointer of the received properties for each near processor + void * ptr_prp = (unsigned char *) recv_mem_gm.get(i).getPointer() + n_ele * sizeof(Point<dim, St> ); + + PtrMemory * ptr1 = new PtrMemory(ptr_pos, n_ele * sizeof(Point<dim, St> )); + PtrMemory * ptr2 = new PtrMemory(ptr_prp, n_ele * sizeof(prop)); + + // create vector representation to a piece of memory already allocated + + openfpm::vector<Point<dim, St>, PtrMemory, typename memory_traits_lin<Point<dim, St>>::type, memory_traits_lin ,openfpm::grow_policy_identity> vpos; + openfpm::vector<prp_object, PtrMemory, typename memory_traits_lin<prp_object>::type, memory_traits_lin ,openfpm::grow_policy_identity> vprp; + + vpos.setMemory(*ptr1); + vprp.setMemory(*ptr2); + + vpos.resize(n_ele); + vprp.resize(n_ele); + + // Add the received particles to v_pos and v_prp + + // source object type + typedef encapc<1, prop, typename openfpm::vector<prop>::layout_type> encap_src; + // destination object type + typedef encapc<1, prp_object, typename openfpm::vector<prp_object>::layout_type> encap_dst; + + size_t j = 0; + for (; j < vpos.size() && o_p_id < out_part.size(); j++, o_p_id++) + { + v_pos.set(out_part.get(o_p_id), vpos.get(j)); + v_prp.template set_o<decltype(vprp.get(j)), prp... >(out_part.get(o_p_id), vprp.get(j)); + } + + for (; j < vpos.size(); j++) + { + v_pos.add(); + v_pos.set(v_pos.size() - 1, vpos.get(j)); + v_prp.template set_o<decltype(vprp.get(j)), prp... >(v_prp.size() - 1, vprp.get(j)); + } + } + + // remove the (out-going particles) in the vector + + v_pos.remove(out_part, o_p_id); + v_prp.remove(out_part, o_p_id); + } + /*! \brief Calculate send buffers total size and allocation * * \tparam prp_object object containing only the properties to send @@ -832,6 +960,86 @@ public: g_m = v_pos.size(); } + /*! \brief It move all the particles that does not belong to the local processor to the respective processor + * + * \tparam out of bound policy it specify what to do when the particles are detected out of bound + * + * In general this function is called after moving the particles to move the + * elements out the local processor. Or just after initialization if each processor + * contain non local particles + * + */ + template<typename obp = KillParticle,unsigned int ... prp> void map_list() + { + // outgoing particles-id + openfpm::vector<size_t> out_part; + + // 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); + + // Contain the processor id of each particle (basically where they have to go) + labelParticleProcessor<obp>(opart, prc_sz, out_part); + + // Calculate the sending buffer size for each processor, put this information in + // a contiguous buffer + p_map_req.resize(v_cl.getProcessingUnits()); + openfpm::vector<size_t> prc_sz_r; + openfpm::vector<size_t> prc_r; + + for (size_t i = 0; i < v_cl.getProcessingUnits(); i++) + { + if (prc_sz.get(i) != 0) + { + p_map_req.get(i) = prc_r.size(); + prc_r.add(i); + prc_sz_r.add(prc_sz.get(i)); + } + } + + // Sending property object + typedef object<typename object_creator<typename prop::type, prp...>::type> prp_object; + + // Allocate the send buffers + + openfpm::vector<pos_prop_sel<prp_object>> pb; + + // fill the send buffers + fill_send_map_buf_list<prp_object,prp...>(prc_r, prc_sz_r, pb); + + // Create the set of pointers + openfpm::vector<void *> ptr(prc_r.size()); + for (size_t i = 0; i < prc_r.size(); i++) + { + ptr.get(i) = pb.get(i).pos.getPointer(); + } + + // convert the particle number to buffer size + for (size_t i = 0; i < prc_sz_r.size(); i++) + { + prc_sz_r.get(i) = prc_sz_r.get(i) * (sizeof(prop) + sizeof(Point<dim, St> )); + } + + // Send and receive the particles + + recv_mem_gm.clear(); + v_cl.sendrecvMultipleMessagesNBX(prc_sz_r.size(), (size_t *) prc_sz_r.getPointer(), (size_t *) prc_r.getPointer(), (void **) ptr.getPointer(), vector_dist::message_alloc_map, this, NEED_ALL_SIZE); + + // Process the incoming particles + + process_received_map_list<prp_object, prp...>(out_part); + + // mark the ghost part + + g_m = v_pos.size(); + } + /*! \brief It synchronize the properties and position of the ghost particles * * \tparam prp list of properties to get synchronize diff --git a/src/Vector/vector_dist_unit_test.hpp b/src/Vector/vector_dist_unit_test.hpp index 35259073aaff8f47401ff58886ab2f0022fcaf7e..68921c245a3281384632d925a0920bf577c1495f 100644 --- a/src/Vector/vector_dist_unit_test.hpp +++ b/src/Vector/vector_dist_unit_test.hpp @@ -60,7 +60,7 @@ template<unsigned int dim> size_t total_n_part_lc(vector_dist<dim,float, Point_t * \param n_out out of domain + ghost particles counter * */ -template<unsigned int dim> inline void count_local_n_local(vector_dist<dim,float, Point_test<float>, CartDecomposition<dim,float> > & vd, vector_dist_iterator & it, size_t (& bc)[dim] , Box<dim,float> & box, Box<dim,float> & dom_ext, size_t & l_cnt, size_t & nl_cnt, size_t & n_out) +template<unsigned int dim,typename vector_dist> inline void count_local_n_local(vector_dist & vd, vector_dist_iterator & it, size_t (& bc)[dim] , Box<dim,float> & box, Box<dim,float> & dom_ext, size_t & l_cnt, size_t & nl_cnt, size_t & n_out) { typedef Point<dim,float> s; const CartDecomposition<dim,float> & ct = vd.getDecomposition(); @@ -789,6 +789,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_map ) } } + BOOST_AUTO_TEST_CASE( vector_dist_not_periodic_map ) { typedef Point<3,float> s; @@ -1158,6 +1159,111 @@ BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test ) } } +BOOST_AUTO_TEST_CASE( vector_dist_periodic_map_list ) +{ + typedef Point<3,float> s; + + Vcluster & v_cl = create_vcluster(); + + if (v_cl.getProcessingUnits() > 3) + return; + + // 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; + + BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector with map_list k=" << k ); + + //! [Create a vector of random elements on each processor 3D] + + Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0}); + + // Boundary conditions + size_t bc[3]={PERIODIC,PERIODIC,PERIODIC}; + + // factor + float factor = pow(create_vcluster().getProcessingUnits()/2.0f,1.0f/3.0f); + + // ghost + Ghost<3,float> ghost(0.05 / factor); + + // ghost2 (a little bigger because of round off error) + Ghost<3,float> ghost2(0.05001 / factor); + + // Distributed vector + vector_dist<3,float, aggregate<float,float,std::list<int>,openfpm::vector<size_t>>, CartDecomposition<3,float> > vd(k,box,bc,ghost); + + 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); + + ++it; + } + + vd.map_list<KillParticle,0,1>(); + + // sync the ghost + vd.ghost_get<0>(); + + //! [Create a vector of random elements on each processor 3D] + + // Domain + ghost + Box<3,float> dom_ext = box; + dom_ext.enlarge(ghost2); + + // Iterate on all particles domain + ghost + size_t l_cnt = 0; + size_t nl_cnt = 0; + size_t n_out = 0; + + auto it2 = vd.getIterator(); + count_local_n_local(vd,it2,bc,box,dom_ext,l_cnt,nl_cnt,n_out); + + // No particles should be out of domain + ghost + BOOST_REQUIRE_EQUAL(n_out,0ul); + + // Ghost must populated because we synchronized them + if (k > 524288) + { + BOOST_REQUIRE(nl_cnt != 0); + BOOST_REQUIRE(l_cnt > nl_cnt); + } + + // Sum all the particles inside the domain + v_cl.sum(l_cnt); + v_cl.execute(); + BOOST_REQUIRE_EQUAL(l_cnt,(size_t)k); + + l_cnt = 0; + nl_cnt = 0; + + // Iterate only on the ghost particles + auto itg = vd.getGhostIterator(); + count_local_n_local(vd,itg,bc,box,dom_ext,l_cnt,nl_cnt,n_out); + + // No particle on the ghost must be inside the domain + BOOST_REQUIRE_EQUAL(l_cnt,0ul); + + // Ghost must be populated + if (k > 524288) + { + BOOST_REQUIRE(nl_cnt != 0); + } +} + BOOST_AUTO_TEST_SUITE_END() #endif /* VECTOR_DIST_UNIT_TEST_HPP_ */