Commit 3099e276 authored by incardon's avatar incardon

Vector with secondary grid distribution working

parent cace612b
......@@ -129,6 +129,9 @@ protected:
//! Structure that store the cartesian grid information
grid_sm<dim, void> gr;
//! Structure that store the cartesian grid information
grid_sm<dim, void> gr_dist;
//! Structure that decompose your structure into cell without creating them
//! useful to convert positions to CellId or sub-domain id in this case
CellDecomposer_sm<dim, T, shift<dim,T>> cd;
......@@ -139,6 +142,10 @@ protected:
//! Box Spacing
T spacing[dim];
//! Magnification factor between distribution and
//! decomposition
size_t magn[dim];
//! Runtime virtual cluster machine
Vcluster & v_cl;
......@@ -166,14 +173,15 @@ protected:
/*! \brief It convert the box from the domain decomposition into sub-domain
*
* The decomposition box from the domain-decomposition contain the box in integer
* coordinates
* coordinates. This box is converted into a continuos box. It also adjust loc_box
* if the distribution grid and the decomposition grid are different.
*
* \param loc_box local box
*
* \return the corresponding sib-domain
* \return the corresponding sub-domain
*
*/
SpaceBox<dim,T> convertDecBoxIntoSubDomain(const SpaceBox<dim,size_t> & loc_box)
template<typename Memory_bx> SpaceBox<dim,T> convertDecBoxIntoSubDomain(encapc<1,::Box<dim,size_t>,Memory_bx> loc_box)
{
// A point with all coordinate to one
size_t one[dim];
......@@ -182,6 +190,15 @@ protected:
SpaceBox<dim, size_t> sub_dc = loc_box;
SpaceBox<dim, size_t> sub_dce = sub_dc;
sub_dce.expand(one);
sub_dce.mul(magn);
// shrink by one
for (size_t i = 0 ; i < dim ; i++)
{
loc_box.template get<Box::p1>()[i] = sub_dce.getLow(i);
loc_box.template get<Box::p2>()[i] = sub_dce.getHigh(i) - 1;
}
SpaceBox<dim, T> sub_d(sub_dce);
sub_d.mul(spacing);
sub_d += domain.getP1();
......@@ -195,7 +212,7 @@ protected:
// domain (avoiding rounding off error)
for (size_t i = 0; i < dim; i++)
{
if (sub_dc.getHigh(i) == cd.getGrid().size(i) - 1)
if (sub_dc.getHigh(i) == gr.size(i) - 1)
sub_d.setHigh(i, domain.getHigh(i));
if (sub_dc.getLow(i) == 0)
......@@ -235,7 +252,7 @@ public:
// Optimize the decomposition creating bigger spaces
// And reducing Ghost over-stress
dec_optimizer<dim, Graph_CSR<nm_v, nm_e>> d_o(dist.getGraph(), gr.getSize());
dec_optimizer<dim, Graph_CSR<nm_v, nm_e>> d_o(dist.getGraph(), gr_dist.getSize());
// Ghost
Ghost<dim,long int> ghe;
......@@ -250,15 +267,16 @@ public:
// optimize the decomposition
d_o.template optimize<nm_v::sub_id, nm_v::proc_id>(dist.getGraph(), p_id, loc_box, box_nn_processor,ghe,bc);
// Initialize ss_box and bbox
// Initialize
if (loc_box.size() >= 0)
{
bbox = convertDecBoxIntoSubDomain(loc_box.get(0));
proc_box = loc_box.get(0);
sub_domains.add(bbox);
}
// convert into sub-domain
for (size_t s = 0; s < loc_box.size(); s++)
for (size_t s = 1; s < loc_box.size(); s++)
{
SpaceBox<dim,T> sub_d = convertDecBoxIntoSubDomain(loc_box.get(s));
......@@ -277,7 +295,8 @@ public:
// fine_s structure contain the processor id for each sub-sub-domain
// with sub-sub-domain we mean the sub-domain decomposition before
// running dec_optimizer (before merging sub-domains)
auto it = dist.getGraph().getVertexIterator();
/* auto it = dist.getGraph().getVertexIterator();
while (it.isNext())
{
......@@ -287,6 +306,24 @@ public:
fine_s.get(key) = dist.getGraph().template vertex_p<nm_v::proc_id>(key);
++it;
}*/
grid_key_dx_iterator<dim> git(gr);
while (git.isNext())
{
auto key = git.get();
grid_key_dx<dim> key2;
for (size_t i = 0 ; i < dim ; i++)
key2.set_d(i,key.get(i) / magn[i]);
size_t lin = gr_dist.LinId(key2);
size_t lin2 = gr.LinId(key);
fine_s.get(lin2) = dist.getGraph().template vertex_p<nm_v::proc_id>(lin);
++git;
}
Initialize_geo_cell_lists();
......@@ -954,15 +991,43 @@ public:
return bc;
}
/*! \brief Calculate magnification
*
* \param gm distribution grid
*
*/
void calculate_magn(const grid_sm<dim,void> & gm)
{
if (gm.size() == 0)
{
for (size_t i = 0 ; i < dim ; i++)
magn[i] = 1;
}
else
{
for (size_t i = 0 ; i < dim ; i++)
{
if (gr.size(i) % gm.size(i) != 0)
std::cerr << __FILE__ << ":" << __LINE__ << ".Error the distribution grid must be multiple of the decomposition grid" << std::endl;
magn[i] = gr.size(i) / gm.size(i);
}
}
}
/*! \brief Set the parameter of the decomposition
*
* \param div_ storing into how many sub-sub-domains to decompose on each dimension
* \param domain_ domain to decompose
* \param bc boundary conditions
* \param ghost Ghost size
* \param sec_dist Distribution grid. The distribution grid help in reducing the underlying
* distribution problem simplifying decomposition problem. This is done in order to
* reduce the load/balancing dynamic load balancing problem
*
*/
void setParameters(const size_t (& div_)[dim], ::Box<dim,T> domain_, const size_t (& bc)[dim] ,const Ghost<dim,T> & ghost)
void setParameters(const size_t (& div_)[dim], ::Box<dim,T> domain_, const size_t (& bc)[dim] ,const Ghost<dim,T> & ghost, const grid_sm<dim,void> & sec_dist = grid_sm<dim,void>())
{
// set the boundary conditions
for (size_t i = 0 ; i < dim ; i++)
......@@ -976,8 +1041,20 @@ public:
domain = domain_;
cd.setDimensions(domain, div_, 0);
// We we have a secondary grid costruct a reduced graph
if (sec_dist.size(0) != 0)
{
calculate_magn(sec_dist);
gr_dist.setDimensions(sec_dist.getSize());
}
else
{
calculate_magn(sec_dist);
gr_dist = gr;
}
// init distribution
dist.createCartGraph(gr, domain);
dist.createCartGraph(gr_dist, domain);
}
......@@ -1397,7 +1474,7 @@ public:
return ghost;
}
/*! \brief Method to access to the grid information of the decomposition
/*! \brief Decomposition grid
*
* \return the grid
*/
......@@ -1406,6 +1483,15 @@ public:
return gr;
}
/*! \brief Distribution grid
*
* \return the grid
*/
const grid_sm<dim,void> getDistGrid()
{
return gr_dist;
}
////////////// Functions to get decomposition information ///////////////
/*! \brief Write the decomposition as VTK file
......
......@@ -342,6 +342,78 @@ BOOST_AUTO_TEST_CASE( CartDecomposition_ext_non_periodic_test)
}
BOOST_AUTO_TEST_CASE( CartDecomposition_non_periodic_test_dist_grid)
{
// Vcluster
Vcluster & vcl = create_vcluster();
CartDecomposition<3, float> dec(vcl);
// Physical domain
Box<3, float> box( { 0.0, 0.0, 0.0 }, { 1.0, 1.0, 1.0 });
size_t div[3];
size_t div_sub[3];
// Get the number of processor and calculate the number of sub-domain
// for each processor (SUB_UNIT_FACTOR=64)
size_t n_proc = vcl.getProcessingUnits();
size_t n_sub = n_proc * SUB_UNIT_FACTOR*4*4*4;
// Set the number of sub-domains on each dimension (in a scalable way)
for (int i = 0; i < 3; i++)
{div[i] = openfpm::math::round_big_2(pow(n_sub,1.0/3));}
// create a sub_distribution grid
for (int i = 0; i < 3; i++)
{div_sub[i] = div[i] / 4;}
grid_sm<3,void> gsub(div_sub);
// Define ghost
Ghost<3, float> g(0.01);
// Boundary conditions
size_t bc[] = { NON_PERIODIC, NON_PERIODIC, NON_PERIODIC };
// Decompose
dec.setParameters(div,box,bc,g,gsub);
dec.decompose();
dec.write("Test_sub_dist2");
// For each calculated ghost box
for (size_t i = 0; i < dec.getNIGhostBox(); i++)
{
SpaceBox<3,float> b = dec.getIGhostBox(i);
size_t proc = dec.getIGhostBoxProcessor(i);
// sample one point inside the box
Point<3,float> p = b.rnd();
// Check that ghost_processorsID return that processor number
const openfpm::vector<size_t> & pr = dec.ghost_processorID<CartDecomposition<3,float>::processor_id>(p);
bool found = false;
for (size_t j = 0; j < pr.size(); j++)
{
if (pr.get(j) == proc)
{ found = true; break;}
}
if (found == false)
{
const openfpm::vector<size_t> pr2 = dec.ghost_processorID<CartDecomposition<3,float>::processor_id>(p);
}
BOOST_REQUIRE_EQUAL(found,true);
}
// Check the consistency
bool val = dec.check_consistency();
BOOST_REQUIRE_EQUAL(val,true);
}
BOOST_AUTO_TEST_SUITE_END()
#endif
......@@ -340,7 +340,7 @@ public:
* ghost communications.
*
*/
vector_dist(size_t np, Box<dim, St> box, const size_t (&bc)[dim], const Ghost<dim, St> & g, size_t opt = 0)
vector_dist(size_t np, Box<dim, St> box, const size_t (&bc)[dim], const Ghost<dim, St> & g, size_t opt = 0, const grid_sm<dim,void> & gdist = grid_sm<dim,void>())
:v_cl(create_vcluster()),opt(opt) SE_CLASS3_VDIST_CONSTRUCTOR
{
#ifdef SE_CLASS2
......@@ -353,7 +353,7 @@ public:
check_parameters(box);
init_structures(np);
this->init_decomposition(box,bc,g,opt);
this->init_decomposition(box,bc,g,opt,gdist);
#ifdef SE_CLASS3
se3.Initialize();
......
......@@ -345,7 +345,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_multiphase_cell_list_sym_test )
verlet_type NNver_all[4];
// This create a Verlet-list between phase all phases to all the other phases
// This create a Verlet-list between each phase to all the other phases
NNver_all[0] = createVerletSymM<2>(0,phases.get(0),phases,CL_all,r_cut);
NNver_all[1] = createVerletSymM<2>(1,phases.get(1),phases,CL_all,r_cut);
NNver_all[2] = createVerletSymM<2>(2,phases.get(2),phases,CL_all,r_cut);
......
......@@ -1186,63 +1186,13 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list_no_bottom )
}
}
BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list )
template<typename part_prop> void test_crs_full(vector_dist<3,float, part_prop > & vd,
vector_dist<3,float, part_prop > & vd2,
std::default_random_engine & eg,
std::uniform_real_distribution<float> & ud,
size_t start,
float r_cut)
{
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 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);
Ghost<3,float> ghost2(r_cut);
ghost2.setLow(0,0.0);
ghost2.setLow(1,0.0);
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())
......@@ -1259,13 +1209,13 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list )
// Fill some properties randomly
vd.getProp<0>(key) = 0;
vd.getProp<1>(key) = 0;
vd.getProp<2>(key) = key.getKey() + start;
vd.template getProp<0>(key) = 0;
vd.template getProp<1>(key) = 0;
vd.template getProp<2>(key) = key.getKey() + start;
vd2.getProp<0>(key) = 0;
vd2.getProp<1>(key) = 0;
vd2.getProp<2>(key) = key.getKey() + start;
vd2.template getProp<0>(key) = 0;
vd2.template getProp<1>(key) = 0;
vd2.template getProp<2>(key) = key.getKey() + start;
++it;
}
......@@ -1273,9 +1223,11 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list )
vd.map();
vd2.map();
Vcluster & v_cl = create_vcluster();
// sync the ghost
vd.ghost_get<0,2>();
vd2.ghost_get<0,2>();
vd.template ghost_get<0,2>();
vd2.template ghost_get<0,2>();
auto NN = vd.getVerlet(r_cut);
auto p_it = vd.getDomainIterator();
......@@ -1286,6 +1238,13 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list )
Point<3,float> xp = vd.getPos(p);
if (v_cl.getProcessUnitID() == 2 && p.getKey() == 137)
{
int debug = 0;
debug++;
}
auto Np = NN.getNNIterator(p.getKey());
while (Np.isNext())
......@@ -1309,10 +1268,10 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list )
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);
vd.template getProp<0>(p)++;
vd.template getProp<3>(p).add();
vd.template getProp<3>(p).last().xq = xq;
vd.template getProp<3>(p).last().id = vd.template getProp<2>(q);
}
++Np;
......@@ -1334,7 +1293,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list )
Point<3,float> xp = vd2.getPos(p);
auto Np = NN2.getNNIterator<NO_CHECK>(p);
auto Np = NN2.template getNNIterator<NO_CHECK>(p);
while (Np.isNext())
{
......@@ -1355,16 +1314,16 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list )
if (distance < r_cut )
{
vd2.getProp<1>(p)++;
vd2.getProp<1>(q)++;
vd2.template getProp<1>(p)++;
vd2.template getProp<1>(q)++;
vd2.getProp<4>(p).add();
vd2.getProp<4>(q).add();
vd2.template getProp<4>(p).add();
vd2.template 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);
vd2.template getProp<4>(p).last().xq = xq;
vd2.template getProp<4>(q).last().xq = xp;
vd2.template getProp<4>(p).last().id = vd2.template getProp<2>(q);
vd2.template getProp<4>(q).last().id = vd2.template getProp<2>(p);
}
++Np;
......@@ -1373,8 +1332,8 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list )
++p_it2;
}
vd2.ghost_put<add_,1>();
vd2.ghost_put<merge_,4>();
vd2.template ghost_put<add_,1>();
vd2.template ghost_put<merge_,4>();
auto p_it3 = vd.getDomainIterator();
......@@ -1383,15 +1342,21 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list )
{
auto p = p_it3.get();
ret &= vd2.getProp<1>(p) == vd.getProp<0>(p);
ret &= vd2.template getProp<1>(p) == vd.template getProp<0>(p);
vd.getProp<3>(p).sort();
vd2.getProp<4>(p).sort();
if (ret == false)
{
Point<3,float> xp = vd2.getPos(p);
std::cout << "ERROR " << vd2.template getProp<1>(p) << " " << vd.template getProp<0>(p) << " " << xp.toString() << std::endl;
}
ret &= vd.getProp<3>(p).size() == vd2.getProp<4>(p).size();
vd.template getProp<3>(p).sort();
vd2.template getProp<4>(p).sort();
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;
ret &= vd.template getProp<3>(p).size() == vd2.template getProp<4>(p).size();
for (size_t i = 0 ; i < vd.template getProp<3>(p).size() ; i++)
ret &= vd.template getProp<3>(p).get(i).id == vd2.template getProp<4>(p).get(i).id;
if (ret == false)
break;
......@@ -1402,6 +1367,149 @@ BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list )
BOOST_REQUIRE_EQUAL(ret,true);
}
void test_csr_verlet_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::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 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);
Ghost<3,float> ghost2(r_cut);
ghost2.setLow(0,0.0);
ghost2.setLow(1,0.0);
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);
test_crs_full(vd,vd2,eg,ud,start,r_cut);
}
void test_csr_verlet_list_override()
{
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 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);
Ghost<3,float> ghost2(r_cut);
ghost2.setLow(0,0.0);
ghost2.setLow(1,0.0);
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;
size_t gdist_d[3];
size_t gdist2_d[3];
gdist_d[0] = 1;
gdist_d[1] = 2;
gdist_d[2] = 5;
gdist2_d[0] = 1;
gdist2_d[1] = 2;
gdist2_d[2] = 5;
grid_sm<3,void> gdist(gdist_d);
grid_sm<3,void> gdist2(gdist2_d);
// Distributed vector
vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_GHOST,gdist_d);
vector_dist<3,float, part_prop > vd2(k,box,bc,ghost2,BIND_DEC_TO_GHOST,gdist2_d);
size_t start = vd.init_size_accum(k);
test_crs_full(vd,vd2,eg,ud,start,r_cut);
}
BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list )
{
test_csr_verlet_list();
}
BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list_dec_override )
{
test_csr_verlet_list_override();
}
BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list_partit )
{
Vcluster & v_cl = create_vcluster();
......
......@@ -881,7 +881,11 @@ public:
* \param opt additional options
*
*/
void init_decomposition(Box<dim,St> & box, const size_t (& bc)[dim],const Ghost<dim,St> & g, size_t opt)
void init_decomposition(Box<dim,St> & box,
const size_t (& bc)[dim],
const Ghost<dim,St> & g,
size_t opt,
const grid_sm<dim,void> & gdist)
{
size_t div[dim];
......@@ -917,7 +921,7 @@ public:
}
// Create the sub-domains
dec.setParameters(div, box, bc, g);
dec.setParameters(div, box, bc, g, gdist);
dec.decompose();
}
......
......@@ -182,14 +182,9 @@ void Test2D_ghost(Box<2,float> & box)