Commit 9ea8975f authored by incardon's avatar incardon

Symmetric Cell-list on progress

parent 6f3a9c09
......@@ -200,7 +200,6 @@ class grid_dist_id
{
// Get the internal ghost boxes and transform into grid units
::Box<dim,St> ib_dom = dec.getProcessorIGhostBox(i,j);
ib_dom -= cd_sm.getOrig();
::Box<dim,long int> ib = cd_sm.convertDomainSpaceIntoGridUnits(ib_dom,dec.periodicity());
// Check if ib is valid if not it mean that the internal ghost does not contain information so skip it
......@@ -321,7 +320,6 @@ class grid_dist_id
{
// Get the internal ghost boxes and transform into grid units
::Box<dim,St> ib_dom = dec.getLocalIGhostBox(i,j);
ib_dom -= cd_sm.getOrig();
::Box<dim,long int> ib = cd_sm.convertDomainSpaceIntoGridUnits(ib_dom,dec.periodicity());
// Check if ib is valid if not it mean that the internal ghost does not contain information so skip it
......
......@@ -78,7 +78,7 @@ BOOST_AUTO_TEST_CASE( grid_dist_id_domain_grid_unit_converter3D_test)
{
// Get the local hyper-cube
SpaceBox<3,float> sub = dec.getSubDomain(i);
sub -= domain.getP1();
// sub -= domain.getP1();
Box<3,size_t> g_box = g_dist.getCellDecomposer().convertDomainSpaceIntoGridUnits(sub,bc);
......
......@@ -78,9 +78,7 @@ template<int dim, typename Decomposition> inline void create_gdb_ext(openfpm::ve
// Get the local sub-domain (Grid conversion must be done with the domain P1 equivalent to 0.0)
// consider that the sub-domain with point P1 equivalent to the domain P1 is a (0,0,0) in grid unit
SpaceBox<Decomposition::dims, typename Decomposition::stype> sp = dec.getSubDomain(i);
sp -= cd_sm.getOrig();
SpaceBox<Decomposition::dims, typename Decomposition::stype> sp_g = dec.getSubDomainWithGhost(i);
sp_g -= cd_sm.getOrig();
// Convert from SpaceBox<dim,St> to SpaceBox<dim,long int>
SpaceBox<Decomposition::dims,long int> sp_t = cd_sm.convertDomainSpaceIntoGridUnits(sp,dec.periodicity());
......
......@@ -347,6 +347,48 @@ public:
return v_prp.template get<id>(g_m - 1);
}
/*! \brief Construct a cell list symmetric based on a cut of radius
*
* \tparam CellL CellList type to construct
*
* \param r_cut interation radius, or size of each cell
*
* \return the Cell list
*
*/
template<typename CellL = CellList<dim, St, FAST, shift<dim, St> > > CellL getCellListSym(St r_cut)
{
// Cell list
CellL cell_list;
size_t div[dim];
// Calculate the Cell list division for this CellList
CellDecomposer<dim,St,shift<dim,St>> cd_sm;
for (size_t i = 0 ; i < dim ; i++)
div[i] = (domain.getHigh(i) - domain.getLow(i)) / r_cut;
size_t pad = 0;
Ghost g = getDecomposition().getGhost();
g.magnify(1.013);
// Calculate the maximum padding
for (size_t i = 0 ; i < dim ; i++)
{
size_t tmp = std::ceil(fabs(g.getLow(i)) / r_cut);
pad = (pad > tmp)?pad:tmp;
}
cd_sm.Initialize(domain,div,pad);
// get the processor bounding box
Box<dim, St> pbox = getDecomposition().getProcessorBounds();
cell_list.Initialize(cd_sm, pbox);
updateCellList(cell_list);
}
/*! \brief Construct a cell list starting from the stored particles
*
* \tparam CellL CellList type to construct
......
......@@ -1817,6 +1817,140 @@ BOOST_AUTO_TEST_CASE( vector_of_vector_dist )
BOOST_REQUIRE_EQUAL(cnt,4*4096ul);
}
BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
{
Vcluster & v_cl = create_vcluster();
if (v_cl.getProcessingUnits() > 24)
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;
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({0.0,0.0,0.0},{1.0,1.0,1.0});
// Boundary conditions
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
float r_cut = 0.1;
// ghost
Ghost<3,float> ghost(r_cut);
typedef aggregate<size_t,size_t> part_prop;
// Distributed vector
vector_dist<3,float, part_prop > 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);
// Fill some properties randomly
vd.getProp<0>(key) = 0.0;
++it;
}
vd.map();
// sync the ghost
vd.ghost_get<0>();
auto NN = vd.getCellList(0.1);
auto p_it = vd.getDomainIterator();
while (p_it.isNext())
{
auto p = p_it.get();
Point<3,float> xp = vd.getPos(p);
auto Np = NN.getIterator(NN.getCell(vd.getPos(p)));
while (Np.isNext())
{
auto q = Np.get();
// 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)++;
++Np;
}
++p_it;
}
// We now try symmetric Cell-list
auto NN2 = vd.getCellListSym(0.1);
auto p_it2 = vd.getDomainIterator();
while (p_it2.isNext())
{
auto p = p_it2.get();
Point<3,float> xp = vd.getPos(p);
auto Np = NN2.getIterator(NN2.getCell(vd.getPos(p)));
while (Np.isNext())
{
auto q = Np.get();
// 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)++;
}
++Np;
}
++p_it;
}
vd.ghost_put<add,1>();
}
#include "vector_dist_cell_list_tests.hpp"
#include "vector_dist_NN_tests.hpp"
#include "vector_dist_complex_prp_unit_test.hpp"
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment