diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp index ded398c37790ab18cb020533efe4ab45e7b36100..03b717856fa4a458f8d54b368fa19a11131f858b 100644 --- a/src/Vector/vector_dist.hpp +++ b/src/Vector/vector_dist.hpp @@ -210,9 +210,13 @@ public: * \param box domain where the vector of elements live * \param bc boundary conditions * \param g Ghost margins + * \param opt additional options. + * * BIND_DEC_TO_GHOST Bind the decomposition to be multiple of the + * ghost size. This is required if we want to use symmetric to eliminate + * ghost communications. * */ - vector_dist(size_t np, Box<dim, St> box, const size_t (&bc)[dim], const Ghost<dim, St> & g) + vector_dist(size_t np, Box<dim, St> box, const size_t (&bc)[dim], const Ghost<dim, St> & g, size_t opt = 0) :v_cl(create_vcluster()) { #ifdef SE_CLASS2 @@ -220,7 +224,7 @@ public: #endif init_structures(np); - this->init_decomposition(box,bc,g); + this->init_decomposition(box,bc,g,opt); } ~vector_dist() diff --git a/src/Vector/vector_dist_cell_list_tests.hpp b/src/Vector/vector_dist_cell_list_tests.hpp index fd599b48e9b3b2930b9b3cf828c16d5a726d09ca..a7dda06b0fc10cdf81b0f8742632f9c1fdfc4ccf 100644 --- a/src/Vector/vector_dist_cell_list_tests.hpp +++ b/src/Vector/vector_dist_cell_list_tests.hpp @@ -735,5 +735,220 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list ) BOOST_REQUIRE_EQUAL(ret,true); } +BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list_no_bottom ) +{ + Vcluster & v_cl = create_vcluster(); + + if (v_cl.getProcessingUnits() > 24) + return; + + float L = 1000.0; + + // set the seed + // create the random generator engine + std::srand(0); + std::default_random_engine eg; + 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 no bottom k=",k); + BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list no bottom 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); + Ghost<3,float> ghost2(r_cut); + ghost2.setLow(2,0.0); + + // 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); + vector_dist<3,float, part_prop > vd2(k,box,bc,ghost2,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); + + vd2.getPos(key)[0] = vd.getPos(key)[0]; + vd2.getPos(key)[1] = vd.getPos(key)[1]; + vd2.getPos(key)[2] = vd.getPos(key)[2]; + + // Fill some properties randomly + + vd.getProp<0>(key) = 0; + vd.getProp<1>(key) = 0; + vd.getProp<2>(key) = key.getKey() + start; + + vd2.getProp<0>(key) = 0; + vd2.getProp<1>(key) = 0; + vd2.getProp<2>(key) = key.getKey() + start; + + ++it; + } + + vd.map(); + vd2.map(); + + // sync the ghost + vd.ghost_get<0,2>(); + vd2.ghost_get<0,2>(); + + auto NN = vd.getVerlet(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(p.getKey()); + + 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 = vd2.getVerletSym(r_cut); + + auto p_it2 = vd2.getDomainIterator(); + + while (p_it2.isNext()) + { + auto p = p_it2.get(); + + Point<3,float> xp = vd2.getPos(p); + + auto Np = NN2.getNNIterator<NO_CHECK>(p.getKey()); + + while (Np.isNext()) + { + auto q = Np.get(); + + if (p.getKey() == q) + { + ++Np; + continue; + } + + // repulsive + + Point<3,float> xq = vd2.getPos(q); + Point<3,float> f = (xp - xq); + + float distance = f.norm(); + + // Particle should be inside r_cut range + + if (distance < r_cut ) + { + vd2.getProp<1>(p)++; + vd2.getProp<1>(q)++; + + vd2.getProp<4>(p).add(); + vd2.getProp<4>(q).add(); + + vd2.getProp<4>(p).last().xq = xq; + vd2.getProp<4>(q).last().xq = xp; + vd2.getProp<4>(p).last().id = vd2.getProp<2>(q); + vd2.getProp<4>(q).last().id = vd2.getProp<2>(p); + } + + ++Np; + } + + + ++p_it2; + } + + vd2.ghost_put<add_,1>(); + vd2.ghost_put<merge_,4>(); + + auto p_it3 = vd.getDomainIterator(); + + bool ret = true; + while (p_it3.isNext()) + { + auto p = p_it3.get(); + + ret &= vd2.getProp<1>(p) == vd.getProp<0>(p); + + + vd.getProp<3>(p).sort(); + vd2.getProp<4>(p).sort(); + + ret &= vd.getProp<3>(p).size() == vd2.getProp<4>(p).size(); + + for (size_t i = 0 ; i < vd.getProp<3>(p).size() ; i++) + ret &= vd.getProp<3>(p).get(i).id == vd2.getProp<4>(p).get(i).id; + + if (ret == false) + break; + + ++p_it3; + } + + BOOST_REQUIRE_EQUAL(ret,true); +} #endif /* SRC_VECTOR_VECTOR_DIST_CELL_LIST_TESTS_HPP_ */ diff --git a/src/Vector/vector_dist_comm.hpp b/src/Vector/vector_dist_comm.hpp index f6ec2e4b5486f1d3d3d46b94e635e35e96f71e31..0115446da6b5f86301e24098bb0ca053acdb7730 100644 --- a/src/Vector/vector_dist_comm.hpp +++ b/src/Vector/vector_dist_comm.hpp @@ -16,6 +16,8 @@ #define NO_POSITION 1 #define WITH_POSITION 2 +#define BIND_DEC_TO_GHOST 1 + /*! \brief This class is an helper for the communication of vector_dist * * \tparam dim Dimensionality of the space where the elements lives @@ -774,22 +776,42 @@ public: * \param box domain * \param bc boundary conditions * \param g ghost extension + * \param opt additional options * */ - void init_decomposition(Box<dim,St> & box, const size_t (& bc)[dim],const Ghost<dim,St> & g) + void init_decomposition(Box<dim,St> & box, const size_t (& bc)[dim],const Ghost<dim,St> & g, size_t opt) { - // Create a valid decomposition of the space - // Get the number of processor and calculate the number of sub-domain - // for decomposition - size_t n_proc = v_cl.getProcessingUnits(); - size_t n_sub = n_proc * getDefaultNsubsub(); - - // Calculate the maximum number (before merging) of sub-domain on - // each dimension size_t div[dim]; - for (size_t i = 0; i < dim; i++) + + if (opt & BIND_DEC_TO_GHOST) + { + // padding + size_t pad = 0; + + // CellDecomposer + CellDecomposer_sm<dim,St,shift<dim,St>> cd_sm; + + // Calculate the divisions for the symmetric Cell-lists + cl_param_calculateSym<dim,St>(box,cd_sm,g,pad); + + for (size_t i = 0 ; i < dim ; i++) + div[i] = cd_sm.getDiv()[i] - 2*pad; + } + else { - div[i] = openfpm::math::round_big_2(pow(n_sub, 1.0 / dim)); + // Create a valid decomposition of the space + // Get the number of processor and calculate the number of sub-domain + // for decomposition + size_t n_proc = v_cl.getProcessingUnits(); + size_t n_sub = n_proc * getDefaultNsubsub(); + + // Calculate the maximum number (before merging) of sub-domain on + // each dimension + + for (size_t i = 0; i < dim; i++) + { + div[i] = openfpm::math::round_big_2(pow(n_sub, 1.0 / dim)); + } } // Create the sub-domains