Commit 8aaf2009 authored by Pietro Incardona's avatar Pietro Incardona

Fixing periodic boundary conditions

parent fd4345ea
......@@ -710,7 +710,10 @@ public:
// Make the id unique
if (opt == UNIQUE)
{
ids_p.sort();
ids_p.unique();
}
return ids_p;
}
......@@ -756,7 +759,10 @@ public:
// Make the id unique
if (opt == UNIQUE)
ids.unique();
{
ids_p.sort();
ids_p.unique();
}
return ids;
}
......@@ -794,7 +800,10 @@ public:
// Make the id unique
if (opt == UNIQUE)
{
ids_p.sort();
ids_p.unique();
}
return ids_p;
}
......@@ -832,7 +841,10 @@ public:
// Make the id unique
if (opt == UNIQUE)
ids.unique();
{
ids_p.sort();
ids_p.unique();
}
return ids;
}
......
......@@ -108,6 +108,7 @@ class ie_loc_ghost
SpaceBox<dim,T> sub_with_ghost = sub_domains_prc.get(j).bx;
size_t rj = sub_domains_prc.get(j).sub;
// Avoid to intersect the box with itself
if (rj == i && sub_domains_prc.get(j).cmb == zero)
continue;
......
......@@ -85,7 +85,7 @@ class grid_dist_iterator<dim,device_grid,FREE>
size_t g_c;
//! List of the grids we are going to iterate
openfpm::vector<device_grid> & gList;
const openfpm::vector<device_grid> & gList;
//! Extension of each grid: domain and ghost + domain
openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext;
......@@ -134,13 +134,8 @@ class grid_dist_iterator<dim,device_grid,FREE>
* \param gk std::vector of the local grid
*
*/
<<<<<<< HEAD
grid_dist_iterator(openfpm::vector<device_grid> & gk, openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext)
:g_c(0),gList(gk),gdb_ext(gdb_ext),m(0)
=======
grid_dist_iterator(const Vcluster_object_array<device_grid> & gk, openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext, grid_key_dx<dim> stop)
grid_dist_iterator(const openfpm::vector<device_grid> & gk, openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext, grid_key_dx<dim> stop)
:g_c(0),gList(gk),gdb_ext(gdb_ext),stop(stop)
>>>>>>> numerics
{
// Initialize the current iterator
// with the first grid
......@@ -248,11 +243,7 @@ class grid_dist_iterator<dim,device_grid,FIXED>
size_t g_c;
//! List of the grids we are going to iterate
<<<<<<< HEAD
openfpm::vector<device_grid> & gList;
=======
const Vcluster_object_array<device_grid> & gList;
>>>>>>> numerics
const openfpm::vector<device_grid> & gList;
//! Extension of each grid: domain and ghost + domain
const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext;
......@@ -297,11 +288,7 @@ class grid_dist_iterator<dim,device_grid,FIXED>
* \param gk std::vector of the local grid
*
*/
<<<<<<< HEAD
grid_dist_iterator(openfpm::vector<device_grid> & gk, const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext)
=======
grid_dist_iterator(const Vcluster_object_array<device_grid> & gk, const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext)
>>>>>>> numerics
grid_dist_iterator(const openfpm::vector<device_grid> & gk, const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext)
:g_c(0),gList(gk),gdb_ext(gdb_ext)
{
// Initialize the current iterator
......
......@@ -34,7 +34,7 @@ class grid_dist_iterator_sub
size_t g_c;
//! List of the grids we are going to iterate
openfpm::vector<device_grid> & gList;
const openfpm::vector<device_grid> & gList;
//! Extension of each grid: domain and ghost + domain
const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext;
......@@ -142,7 +142,7 @@ class grid_dist_iterator_sub
* \param gdb_ext information about the local grids
*
*/
grid_dist_iterator_sub(const grid_key_dx<dim> & start, const grid_key_dx<dim> & stop ,const Vcluster_object_array<device_grid> & gk, const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext)
grid_dist_iterator_sub(const grid_key_dx<dim> & start, const grid_key_dx<dim> & stop ,const openfpm::vector<device_grid> & gk, const openfpm::vector<GBoxes<device_grid::dims>> & gdb_ext)
:g_c(0),gList(gk),gdb_ext(gdb_ext),start(start),stop(stop)
{
// Initialize the current iterator
......
......@@ -1132,7 +1132,7 @@ void Test3D_decit(const Box<3,float> & domain, long int k)
// Grid sm
grid_sm<3,void> info(sz);
auto dom = g_dist.getSubDomainIterator({0,0,0},{sz[0]-2,sz[1]-2,sz[2]-2});
auto dom = g_dist.getSubDomainIterator({0,0,0},{(long int)sz[0]-2,(long int)sz[1]-2,(long int)sz[2]-2});
bool match = true;
......
......@@ -232,7 +232,6 @@ private:
// add particle to communicate
ghost_prc_sz.get(p_id)++;
opart.get(p_id).add(key);
oshift.get(p_id).add(vp_id.get(i).second);
}
......@@ -289,10 +288,15 @@ private:
// get the shift vectors
const openfpm::vector<Point<dim,St>> & shifts = dec.getShiftVectors();
// this map is used to check if a combination is already present
std::unordered_map<size_t, size_t> map_cmb;
// Add local particles coming from periodic boundary, the only boxes that count are the one
// at the border, filter them
// touching the border, filter them
openfpm::vector_std<Box<dim,St>> box_f;
// The boxes touching the border of the domain are divided in groups (first vector)
// each group contain internal ghost coming from sub-domain of the same section
openfpm::vector_std<openfpm::vector_std<Box<dim,St>>> box_f;
openfpm::vector_std<comb<dim>> box_cmb;
for (size_t i = 0 ; i < dec.getNLocalSub() ; i++)
......@@ -301,12 +305,28 @@ private:
for (size_t j = 0 ; j < Nl ; j++)
{
// If they are not in the border the combination is all zero
// If the ghost does not come from the intersection with an out of
// border sub-domain the combination is all zero and n_zero return dim
if (dec.getLocalIGhostPos(i,j).n_zero() == dim)
continue;
box_f.add(dec.getLocalIGhostBox(i,j));
box_cmb.add(dec.getLocalIGhostPos(i,j));
// Check if we already have boxes with such combination
auto it = map_cmb.find(dec.getLocalIGhostPos(i,j).lin());
if (it == map_cmb.end())
{
// we do not have it
box_f.add();
box_f.last().add(dec.getLocalIGhostBox(i,j));
box_cmb.add(dec.getLocalIGhostPos(i,j));
map_cmb[dec.getLocalIGhostPos(i,j).lin()] = box_f.size()-1;
}
else
{
// we have it
box_f.get(it->second).add(dec.getLocalIGhostBox(i,j));
}
}
}
......@@ -316,7 +336,7 @@ private:
else
{
// Label the internal (assigned) particles
auto it = v_pos.getIterator();
auto it = v_pos.getIteratorTo(g_m);
while (it.isNext())
{
......@@ -325,16 +345,29 @@ private:
// If particles are inside these boxes
for (size_t i = 0 ; i < box_f.size() ; i++)
{
if (box_f.get(i).isInside(v_pos.get(key)) == true)
for (size_t j = 0 ; j < box_f.get(i).size() ; j++)
{
Point<dim,St> p = v_pos.get(key);
// shift
p -= shifts.get(box_cmb.get(i).lin());
// add this particle shifting its position
v_pos.add(p);
v_prp.add();
v_prp.last() = v_prp.get(key);
if (box_f.get(i).get(j).isInside(v_pos.get(key)) == true)
{
Point<dim,St> p = v_pos.get(key);
// shift
p -= shifts.get(box_cmb.get(i).lin());
// add this particle shifting its position
v_pos.add(p);
v_prp.add();
v_prp.last() = v_prp.get(key);
// boxes in one group can be overlapping
// we do not have to search for the other
// boxes otherwise we will have duplicate particles
//
// A small note overlap of boxes across groups is fine
// (and needed) because each group has different shift
// producing non overlapping particles
//
break;
}
}
}
......@@ -995,7 +1028,6 @@ public:
*/
template<typename CellL=CellList<dim,St,FAST,shift<dim,St> > > CellL getCellList(St r_cut, const Ghost<dim,St> & enlarge)
{
CellL cell_list;
// calculate the parameters of the cell list
......@@ -1037,6 +1069,63 @@ public:
return cell_list;
}
/*! \brief for each particle get the verlet list
*
* \param verlet output verlet list for each particle
* \param r_cut cut-off radius
*
*/
void getVerlet(openfpm::vector<openfpm::vector<size_t>> & verlet, St r_cut)
{
// resize verlet to store the number of particles
verlet.resize(size_local());
// get the cell-list
auto cl = getCellList(r_cut);
// square of the cutting radius
St r_cut2 = r_cut*r_cut;
// iterate the particles
auto it_p = this->getDomainIterator();
while (it_p.isNext())
{
// key
vect_dist_key_dx key = it_p.get();
// Get the position of the particles
Point<dim,St> p = this->template getPos<0>(key);
// Clear the neighborhood of the particle
verlet.get(key.getKey()).clear();
// Get the neighborhood of the particle
auto NN = cl.template getNNIterator<NO_CHECK>(cl.getCell(p));
while(NN.isNext())
{
auto nnp = NN.get();
// p != q
if (nnp == key.getKey())
{
++NN;
continue;
}
Point<dim,St> q = this->template getPos<0>(nnp);
if (p.distance2(q) < r_cut2)
verlet.get(key.getKey()).add(nnp);
// Next particle
++NN;
}
// next particle
++it_p;
}
}
/*! \brief It return the number of particles contained by the previous processors
*
* \Warning It only work with the initial decomposition
......@@ -1098,8 +1187,8 @@ public:
start.set_d(i,0);
if (dec.isPeriodic(i) == PERIODIC)
{
sz_g[i] = sz[i]-1;
stop.set_d(i,sz_g[i]-1);
sz_g[i] = sz[i];
stop.set_d(i,sz_g[i]-2);
}
else
{
......
......@@ -92,6 +92,12 @@ template<unsigned int dim> inline void count_local_n_local(vector_dist<dim,float
BOOST_AUTO_TEST_SUITE( vector_dist_test )
void print_test(std::string test, size_t sz)
{
if (global_v_cluster->getProcessUnitID() == 0)
std::cout << test << " " << sz << "\n";
}
BOOST_AUTO_TEST_CASE( vector_dist_ghost )
{
// Communication object
......@@ -1031,100 +1037,138 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_interacting_particles )
BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
{
typedef Point<3,float> s;
size_t k = 64*64*64*global_v_cluster->getProcessingUnits();
k = std::pow(k, 1/3.);
// we create a 128x128x128 Grid iterator
size_t sz[3] = {128,128,128};
size_t total = sz[0]*sz[1]*sz[2];
long int big_step = k / 30;
big_step = (big_step == 0)?1:big_step;
long int small_step = 21;
Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
print_test( "Vector cell and verlet list test k<=",k);
// Boundary conditions
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
// 3D test
for ( ; k >= 2 ; k-= (k > 2*big_step)?big_step:small_step )
{
typedef Point<3,float> s;
// factor
float factor = pow(global_v_cluster->getProcessingUnits()/2.0f,1.0f/3.0f);
print_test( "Vector cell and verlet list test k<=",k);
// ghost
Ghost<3,float> ghost(0.05 / factor);
Vcluster & v_cl = *global_v_cluster;
// Distributed vector
vector_dist<3,float, Point_test<float>, CartDecomposition<3,float> > vd(total,box,bc,ghost);
const size_t Ng = k;
// Put particles on a grid creating a Grid iterator
auto it = vd.getGridIterator(sz);
auto it_p = vd.getDomainIterator();
// we create a 128x128x128 Grid iterator
size_t sz[3] = {Ng,Ng,Ng};
while (it_p.isNext())
{
auto key_p = it_p.get();
auto key = it.get();
Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
vd.template getPos<s::x>(key_p)[0] = key.get(0) * it.getSpacing(0);
vd.template getPos<s::x>(key_p)[1] = key.get(1) * it.getSpacing(1);
vd.template getPos<s::x>(key_p)[2] = key.get(2) * it.getSpacing(2);
// Boundary conditions
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
++it;
++it_p;
}
// ghost
Ghost<3,float> ghost(1.0/(Ng-2));
vd.map();
// Distributed vector
vector_dist<3,float, Point_test<float>, CartDecomposition<3,float> > vd(0,box,bc,ghost);
// calculate the distance of the first, second and third neighborhood particle
// Consider that they are on a regular grid
// Put particles on a grid creating a Grid iterator
auto it = vd.getGridIterator(sz);
float spacing = it.getSpacing(0);
float first_dist = spacing;
float second_dist = sqrt(2.0*spacing*spacing);
float third_dist = sqrt(3.0 * spacing*spacing);
while (it.isNext())
{
vd.add();
// add a 5% to dist
auto key = it.get();
first_dist += first_dist * 0.05;
second_dist += second_dist * 0.05;
third_dist += third_dist * 0.05;
vd.template getLastPos<s::x>()[0] = key.get(0) * it.getSpacing(0);
vd.template getLastPos<s::x>()[1] = key.get(1) * it.getSpacing(1);
vd.template getLastPos<s::x>()[2] = key.get(2) * it.getSpacing(2);
// Create a verlet list for each particle
++it;
}
openfpm::vector<openfpm::vector<size_t>> verlet;
vd.getVerlet(verlet);
// distribute particles and sync ghost
vd.map();
bool correct = true;
// Check that the sum of all the particles is the grid size
size_t total = vd.size_local();
v_cl.sum(total);
v_cl.execute();
// for each particle
for (size_t i = 0 ; i < verlet.size() ; i++)
{
BOOST_REQUIRE_EQUAL(total,(Ng-1) * (Ng-1) * (Ng-1));
vd.ghost_get<0>();
// calculate the distance of the first, second and third neighborhood particle
// Consider that they are on a regular grid
float spacing = it.getSpacing(0);
float first_dist = spacing;
float second_dist = sqrt(2.0*spacing*spacing);
float third_dist = sqrt(3.0 * spacing*spacing);
// add a 5% to dist
first_dist += first_dist * 0.05;
second_dist += second_dist * 0.05;
third_dist += third_dist * 0.05;
// Create a verlet list for each particle
openfpm::vector<openfpm::vector<size_t>> verlet;
// first NN
size_t first_NN = 0;
size_t second_NN = 0;
size_t third_NN = 0;
vd.getVerlet(verlet,third_dist);
Point<3,float> p = vd.getPos<0>(i);
//
if (v_cl.getProcessUnitID() == 0)
{
Point<3,float> p1 = vd.template getPos<0>(verlet.get(0).get(4));
Point<3,float> p2 = vd.template getPos<0>(verlet.get(0).get(5));
std::cout << "Point 1: " << p1.toString() << "\n Point 2: " << p2.toString() << "\n";
}
bool correct = true;
// for each neighborhood particle
for (size_t j = 0 ; j < verlet.get(i).size() ; j++)
// for each particle
for (size_t i = 0 ; i < verlet.size() ; i++)
{
auto & NN = verlet.get(i);
// first NN
size_t first_NN = 0;
size_t second_NN = 0;
size_t third_NN = 0;
Point<3,float> q = vc.getPos<0>(NN.get(j));
Point<3,float> p = vd.getPos<0>(i);
float dist = p.distance(q);
// for each neighborhood particle
for (size_t j = 0 ; j < verlet.get(i).size() ; j++)
{
auto & NN = verlet.get(i);
if (dist <= first_dist)
first_NN++;
else if (dist <= second_dist)
second_NN++;
else
third_NN++;
Point<3,float> q = vd.getPos<0>(NN.get(j));
float dist = p.distance(q);
if (dist <= first_dist)
first_NN++;
else if (dist <= second_dist)
second_NN++;
else
third_NN++;
}
correct &= (first_NN == 6);
correct &= (second_NN == 12);
correct &= (third_NN == 8);
if (correct == false)
{
int debug = 0;
debug++;
}
}
correct &= (first_NN == 6);
correct &= (second_NN == 12);
correct &= (third_NN = 8);
BOOST_REQUIRE_EQUAL(correct,true);
}
BOOST_REQUIRE_EQUAL(correct,true);
}
BOOST_AUTO_TEST_SUITE_END()
......
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