Skip to content

GitLab

O
openfpm_pdata
Commit ac799346 authored by Pietro Incardona's avatar Pietro Incardona
Browse files
Options

Fixing periodic

parent f59eade6
Hide whitespace changes
Inline Side-by-side
Showing with 559 additions and 50 deletions
CHANGELOG.md
View file @ ac799346
......@@ -6,9 +6,11 @@ All notable changes to this project will be documented in this file.
### Added
- Grid with periodic boundary conditions
- VTK Writer for distributed grid, now is the default writer
- Installation of linear algebra packages
### Fixed
- GPU compilation
- PARMetis automated installation
### Changed
......
src/Decomposition/CartDecomposition.hpp
View file @ ac799346
......@@ -70,7 +70,7 @@
* * Near processor sub-domain: is a sub-domain that live in the a near (or contiguous) processor
* * Near processor list: the list of all the near processor of the local processor (each processor has a list
* of the near processor)
* * Local ghosts interal or external are all the ghosts that does not involve inter-processor communications
* * Local ghosts internal or external are all the ghosts that does not involve inter-processor communications
*
* \see calculateGhostBoxes() for a visualization of internal and external ghost boxes
*
......
src/Decomposition/Distribution/ParMetisDistribution.hpp
View file @ ac799346
......@@ -350,8 +350,10 @@ public:
}
}
v_cl.sendrecvMultipleMessagesNBX(prc.size(), &sz.get(0), &prc.get(0), &ptr.get(0), message_receive, &partitions,
NONE);
if (prc.size() == 0)
v_cl.sendrecvMultipleMessagesNBX(0, NULL, NULL, NULL, message_receive, &partitions,NONE);
else
v_cl.sendrecvMultipleMessagesNBX(prc.size(), &sz.get(0), &prc.get(0), &ptr.get(0), message_receive, &partitions,NONE);
// Update graphs with the received data
updateGraphs();
......
src/Decomposition/common.hpp
View file @ ac799346
......@@ -42,12 +42,6 @@ struct Box_loc_sub
:bx(bx),sub(sub),cmb(cmb)
{};
template <typename Memory> Box_loc_sub(const Box_loc_sub<dim,T> & bls)
{
bx = bls.bx;
this->sub = bls.sub;
};
Box_loc_sub operator=(const Box<dim,T> & box)
{
::Box<dim,T>::operator=(box);
......@@ -59,7 +53,7 @@ struct Box_loc_sub
};
/*! It contain a box definition and from witch sub-domain it come from (in the local processor)
* and an unique across adjacent processors (for communication)
* and an unique if across adjacent processors (for communication)
*
* If the box come from the intersection of an expanded sub-domain and a sub-domain
*
......@@ -85,22 +79,17 @@ struct Box_loc_sub
template<unsigned int dim, typename T>
struct Box_sub
{
//! Internal ghost box definition
Box<dim,T> bx;
// Domain id
//! Domain id
size_t sub;
// Id
//! see ebx_ibx_form in ie_ghost for the meaning
size_t id;
Box_sub operator=(const Box<dim,T> & box)
{
bx = box;
return *this;
}
//! see ie_ghost follow sector explanation
comb<dim> cmb;
};
//! Particular case for local internal ghost boxes
......@@ -124,13 +113,6 @@ struct Box_sub_k
cmb.zero();
}
Box_sub_k operator=(const Box<dim,T> & box)
{
bx = box;
return *this;
}
// encap interface to make compatible with OpenFPM_IO
template <int i> auto get() -> decltype( std::declval<Box<dim,T> *>()->template get<i>())
{
......
src/Decomposition/ie_ghost.hpp
View file @ ac799346
......@@ -246,8 +246,9 @@ protected:
for (size_t b = 0 ; b < nn_processor_subdomains_g.size() ; b++)
{
::Box<dim,T> bi;
::Box<dim,T> sub_bb(nn_processor_subdomains_g.get(b));
bool intersect = sub_with_ghost.Intersect(::Box<dim,T>(nn_processor_subdomains_g.get(b)),bi);
bool intersect = sub_with_ghost.Intersect(sub_bb,bi);
if (intersect == true)
{
......@@ -272,8 +273,9 @@ protected:
vb_ext.add(pb);
box_nn_processor_int_gg.add(bi);
proc_int_box_g.ebx.add();
proc_int_box_g.ebx.last() = bi;
proc_int_box_g.ebx.last().bx = bi;
proc_int_box_g.ebx.last().sub = i;
proc_int_box_g.ebx.last().cmb = nnpsg_pos.get(b);
// Search where the sub-domain i is in the sent list for processor p_id
size_t k = link_ebx_ibx(nn_p,p_id,i);
......@@ -374,8 +376,9 @@ protected:
// store the box in proc_int_box storing from which sub-domain they come from
Box_sub<dim,T> sb;
sb = b_int.box;
sb.bx = b_int.box;
sb.sub = i;
sb.cmb = nn_p_box_pos.get(k);
size_t p_idp = nn_p.ProctoID(p_id);
......@@ -411,6 +414,7 @@ protected:
}
}
public:
//! Default constructor
......@@ -555,11 +559,39 @@ public:
return proc_int_box.get(id).ebx.get(j).bx;
}
/*! \brief Get the j External ghost box sector
*
* \param id near processor list id (the id go from 0 to getNNProcessor())
* \param j box (each near processor can produce more than one external ghost box)
* \return the sector
*
*/
inline const comb<dim> & getProcessorEGhostPos(size_t id, size_t j) const
{
return proc_int_box.get(id).ebx.get(j).cmb;
}
/*! \brief Get the ghost box sector of the external ghost box linked with the j internal ghost box
*
* \param id near processor list id (the id go from 0 to getNNProcessor())
* \param j box (each near processor can produce more than one internal ghost box)
* \return the sector
*
*/
inline const comb<dim> & getProcessorIGhostPos(size_t id, size_t j) const
{
return proc_int_box.get(id).ibx.get(j).cmb;
}
/*! \brief Get the j Internal ghost box id
*
* Every internal ghost box has a linked external ghost box, because they overlap
* and they must contain the same information (Think on a ghost_get). So if exist
* an internal ghost box with id x, exist also an external ghost box with id x
*
* \param id near processor list id (the id go from 0 to getNNProcessor())
* \param j box (each near processor can produce more than one internal ghost box)
* \return the box
* \return the box id
*
*/
inline size_t getProcessorIGhostId(size_t id, size_t j) const
......@@ -568,6 +600,10 @@ public:
}
/*! \brief Get the j External ghost box id
*
* Every external ghost box has a linked internal ghost box, because they overlap
* and they must contain the same information (Think on a ghost_get). So if exist
* an internal ghost box with id x, exist also an external ghost box with id x
*
* \param id near processor list id (the id go from 0 to getNNProcessor())
* \param j box (each near processor can produce more than one external ghost box)
......@@ -619,7 +655,7 @@ public:
* \return the internal ghost box
*
*/
inline ::Box<dim,T> getIGhostBox(size_t b_id) const
inline const ::Box<dim,T> & getIGhostBox(size_t b_id) const
{
return vb_int.get(b_id).box;
}
......
src/Decomposition/ie_loc_ghost.hpp
View file @ ac799346
......@@ -68,7 +68,8 @@ class ie_loc_ghost
{
Box_sub<dim,T> b;
b.sub = rj;
b = bi;
b.bx = bi;
b.cmb = sub_domains_prc.get(j).cmb;
// local external ghost box
loc_ghost_box.get(i).ebx.add(b);
......@@ -123,7 +124,7 @@ class ie_loc_ghost
{
Box_sub_k<dim,T> b;
b.sub = rj;
b = bi;
b.bx = bi;
b.k = -1;
b.cmb = sub_domains_prc.get(j).cmb;
......@@ -330,7 +331,7 @@ public:
}
/*! \brief For the sub-domain i intersected with the sub-domain j enlarged, the associated
* external ghost box is located in getLocalIGhostBox(j,k) with
* external ghost box is located in getLocalEGhostBox(j,k) with
* getLocalIGhostSub(j,k) == i, this function return k
*
*
......@@ -427,6 +428,18 @@ public:
return loc_ghost_box.get(i).ebx.get(j).bx;
}
/*! \brief Get the j external local ghost box for the local processor
*
* \param i sub-domain
* \param j box
* \return the box
*
*/
inline const comb<dim> & getLocalEGhostPos(size_t i, size_t j) const
{
return loc_ghost_box.get(i).ebx.get(j).cmb;
}
/*! \brief Considering that sub-domain has N internal local ghost box identified
* with the 0 <= k < N that come from the intersection of 2 sub-domains i and j
* where j is enlarged, given the sub-domain i and the id k, it return the id of
......
src/Decomposition/nn_processor.hpp
View file @ ac799346
......@@ -48,6 +48,81 @@ class nn_prcs
//! applyBC function is suppose to be called only one time
bool aBC;
/*! \brief It shift a box but it does consistently
*
* In calculating internal and external ghost boxes, domains are shifted by periodicity.
* In particular, consider a box touching with the left bolder the left border of the domain
*
before shift after shift
+-----------------------------+ +------------------------------+
| | | |
| domain | | domain |
| | | |
| | | |
+---------+ | | +---------+
| | | | | |
| | | | | |
| box | | | | box |
| | | | | |
| | | | | |
+---------+ | | +---------+
| | | |
| | | |
| | | |
| | | |
| | | |
+-----------------------------+ +------------------------------+
*
*
*
*
*
* shifting the box on the right by the size of the domain, we expect to have a box touching with
* the left side the right side of the domain. Because of rounding off problem this is not possible
* with a simple shift. This function ensure consistency like ensuring the previous condition, with
* the assumption that the shift is +/- the domain size
*
* \param box to shift
* \param domain
* \param shift
*
*/
inline void consistent_shift(Box<dim,T> & box, const Box<dim,T> & domain, const Point<dim,T> & shift)
{
for (size_t k = 0 ; k < dim ; k++)
{
// if it touch on the left and shift on the right
if (box.getLow(k) == domain.getLow(k) && shift.get(k) > 0)
{
box.setLow(k,domain.getHigh(k));
box.setHigh(k,box.getHigh(k) + shift.get(k));
}
else if (box.getLow(k) == domain.getHigh(k) && shift.get(k) < 0)
{
box.setLow(k,domain.getLow(k));
box.setHigh(k,box.getHigh(k) + shift.get(k));
}
else if (box.getHigh(k) == domain.getHigh(k) && shift.get(k) < 0)
{
box.setHigh(k,domain.getLow(k));
box.setLow(k,box.getLow(k) + shift.get(k));
}
else if (box.getHigh(k) == domain.getLow(k) && shift.get(k) > 0)
{
box.setHigh(k,domain.getHigh(k));
box.setLow(k,box.getLow(k) + shift.get(k));
}
else
{
box.setHigh(k,box.getHigh(k) + shift.get(k));
box.setLow(k,box.getLow(k) + shift.get(k));
}
}
}
/*! \brief Message allocation
*
* \param message size required to receive from i
......@@ -151,7 +226,14 @@ class nn_prcs
if (sub.Intersect(bp,b_int) == true)
{
sub += shift;
Box<dim,T> sub2 = sub;
sub2 += shift;
// Here we have to be careful of rounding off problems, in particular if any part
// of the sub-domain touch the border of the domain
consistent_shift(sub,domain,shift);
add_nn_subdomain(IDtoProc(k),l,sub,cmbs[j]);
}
}
......
src/Grid/grid_dist_id.hpp
View file @ ac799346
......@@ -19,6 +19,14 @@
#include "Packer_Unpacker/Unpacker.hpp"
#include "Decomposition/CartDecomposition.hpp"
// External ghost box to send for internal ghost box fixation
template<unsigned int dim>
struct Box_fix
{
Box<dim,size_t> bx;
size_t g_id;
};
#define GRID_SUB_UNIT_FACTOR 64
/*! \brief This is a distributed grid
......@@ -75,8 +83,13 @@ class grid_dist_id
Vcluster & v_cl;
//! It map a global ghost id (g_id) to the external ghost box information
//! It is unique across all the near processor
std::unordered_map<size_t,size_t> g_id_to_external_ghost_box;
//! It map a global ghost id (g_id) to the internal ghost box information
//! (is unique for processor), it is not unique across all the near processor
openfpm::vector<std::unordered_map<size_t,size_t>> g_id_to_internal_ghost_box;
// Receiving size
openfpm::vector<size_t> recv_sz;
......@@ -122,10 +135,12 @@ class grid_dist_id
*/
void create_ig_box()
{
if (init_i_g_box == true) return;
// Get the grid info
auto g = cd_sm.getGrid();
if (init_i_g_box == true) return;
g_id_to_internal_ghost_box.resize(dec.getNNProcessors());
// Get the number of near processors
for (size_t i = 0 ; i < dec.getNNProcessors() ; i++)
......@@ -152,7 +167,10 @@ class grid_dist_id
bid_t.box = cvt;
bid_t.g_id = dec.getProcessorIGhostId(i,j);
bid_t.sub = dec.getProcessorIGhostSub(i,j);
bid_t.cmb = dec.getProcessorIGhostPos(i,j);
pib.bid.add(bid_t);
g_id_to_internal_ghost_box.get(i)[bid_t.g_id] = pib.bid.size()-1;
}
}
......@@ -198,6 +216,8 @@ class grid_dist_id
bid_t.sub = sub_id;
bid_t.g_e_box = ib;
bid_t.l_e_box = ib;
bid_t.cmb = dec.getProcessorEGhostPos(i,j);
bid_t.g_id = dec.getProcessorEGhostId(i,j);
// Translate in local coordinate
Box<dim,long int> tb = ib;
tb -= gdb_ext.get(sub_id).origin;
......@@ -285,6 +305,8 @@ class grid_dist_id
pib.bid.add();
pib.bid.last().box = ib;
pib.bid.last().sub = dec.getLocalEGhostSub(i,j);
pib.bid.last().cmb = dec.getLocalEGhostPos(i,j);
pib.bid.last().cmb.sign_flip();
}
}
......@@ -320,10 +342,14 @@ class grid_dist_id
bx_dst -= gdb_ext.get(sub_id_dst).origin;
// create 2 sub grid iterator
if (bx_dst.isValid() == false)
continue;
grid_key_dx_iterator_sub<dim> sub_src(loc_grid.get(i).getGrid(),bx_src.getKP1(),bx_src.getKP2());
grid_key_dx_iterator_sub<dim> sub_dst(loc_grid.get(sub_id_dst).getGrid(),bx_dst.getKP1(),bx_dst.getKP2());
#ifdef DEBUG
#ifdef SE_CLASS1
if (loc_eg_box.get(sub_id_dst).bid.get(k).sub != i)
std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " source and destination are not correctly linked" << "\n";
......@@ -1087,6 +1113,9 @@ public:
//! id
size_t g_id;
//! Sector where it live the linked external ghost box
comb<dim> cmb;
//! sub
size_t sub;
};
......@@ -1119,11 +1148,17 @@ public:
//! Box defining the external ghost box in local coordinates
::Box<dim,size_t> l_e_box;
//! Sector position of the external ghost
comb<dim> cmb;
//! Id
size_t g_id;
//! sub_id in which sub-domain this box live
size_t sub;
};
/*! \brief It store the information about the external ghost box
/*! \brief It store the information about the local external ghost box
*
*
*/
......@@ -1132,6 +1167,9 @@ public:
//! Box defining the external ghost box in local coordinates
::Box<dim,size_t> box;
//! Sector position of the local external ghost box
comb<dim> cmb;
//! sub_id in which sub-domain this box live
size_t sub;
};
......@@ -1190,6 +1228,9 @@ public:
//! Flag that indicate if the internal ghost box has been initialized
bool init_i_g_box = false;
//! Flag that indicate if the internal and external ghost box has been fixed
bool init_fix_ie_g_box = false;
//! Internal ghost boxes in grid units
openfpm::vector<ip_box_grid> ig_box;
......@@ -1241,6 +1282,10 @@ public:
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;
if (g_ig_box.isValid() == false)
continue;
g_ig_box -= gdb_ext.get(sub_id).origin.template convertPoint<size_t>();
// Pack a size_t for the internal ghost id
......@@ -1270,6 +1315,10 @@ public:
// for each ghost box
for (size_t j = 0 ; j < ig_box.get(i).bid.size() ; j++)
{
// we pack only if it is valid
if (ig_box.get(i).bid.get(j).box.isValid() == false)
continue;
// And linked sub-domain
size_t sub_id = ig_box.get(i).bid.get(j).sub;
// Internal ghost box
......
src/Grid/grid_dist_id_unit_test.hpp
View file @ ac799346
......@@ -673,7 +673,7 @@ void Test3D_gg(const Box<3,float> & domain, long int k, long int gk)
*
*/
void Test3D_domain(const Box<3,float> & domain, long int k)
void Test3D_domain(const Box<3,float> & domain, long int k, const periodicity<3> & pr)
{
long int big_step = k / 30;
big_step = (big_step == 0)?1:big_step;
......@@ -699,7 +699,9 @@ void Test3D_domain(const Box<3,float> & domain, long int k)
Ghost<3,float> g(0.01 / factor);
// Distributed grid with id decomposition
grid_dist_id<3, float, aggregate<long int,long int>, CartDecomposition<3,float>> g_dist(sz,domain,g);
grid_dist_id<3, float, aggregate<long int,long int>, CartDecomposition<3,float>> g_dist(sz,domain,g,pr);
auto & v_cl = create_vcluster();
// check the consistency of the decomposition
bool val = g_dist.getDecomposition().check_consistency();
......@@ -732,7 +734,6 @@ void Test3D_domain(const Box<3,float> & domain, long int k)
// Get the total size of the local grids on each processors
// and the total size
Vcluster & v_cl = g_dist.getVC();
v_cl.sum(count2);
v_cl.allGather(count,pnt);
v_cl.execute();
......@@ -1386,7 +1387,7 @@ void Test3D_periodic(const Box<3,float> & domain, long int k)
periodicity<3> pr = {{PERIODIC,PERIODIC,PERIODIC}};
// Distributed grid with id decomposition
grid_dist_id<3, float, Point_test<float>, CartDecomposition<3,float>> g_dist(sz,domain,g,pr);
grid_dist_id<3, float, aggregate<long int>, CartDecomposition<3,float>> g_dist(sz,domain,g,pr);
// check the consistency of the decomposition
bool val = g_dist.getDecomposition().check_consistency();
......@@ -1397,6 +1398,19 @@ void Test3D_periodic(const Box<3,float> & domain, long int k)
size_t count = 0;
// Set to zero the full grid
auto dom1 = g_dist.getDomainGhostIterator();
while (dom1.isNext())
{
auto key = dom1.get();
g_dist.template get<0>(key) = -1;
++dom1;
}
auto dom = g_dist.getDomainIterator();
while (dom.isNext())
......@@ -1430,8 +1444,12 @@ void Test3D_periodic(const Box<3,float> & domain, long int k)
// Domain + Ghost iterator
auto dom_gi = g_dist.getDomainGhostIterator();
size_t out_cnt = 0;
while (dom_gi.isNext())
{
bool out_p = false;
auto key = dom_gi.get();
auto key_g = g_dist.getGKey(key);
......@@ -1439,18 +1457,27 @@ void Test3D_periodic(const Box<3,float> & domain, long int k)
for (size_t i = 0 ; i < 3 ; i++)
{
if (key_g.get(i) < 0)
key_g.set_d(i,key_g.get(i) + k);
{key_g.set_d(i,key_g.get(i) + k);out_p = true;}
else if (key_g.get(i) >= k)
key_g.set_d(i,key_g.get(i) - k);
{key_g.set_d(i,key_g.get(i) - k);out_p = true;}
}
if ( info.LinId(key_g) != g_dist.template get<0>(key) )
if (g_dist.template get<0>(key) != -1 && out_p == true)
out_cnt++;
if ( g_dist.template get<0>(key) != -1 && info.LinId(key_g) != g_dist.template get<0>(key) )
match &= false;
++dom_gi;
}
BOOST_REQUIRE_EQUAL(match, true);
if (k > 83)
{
vcl.sum(out_cnt);
vcl.execute();
BOOST_REQUIRE(out_cnt != 0ul);
}
}
}
......@@ -1544,9 +1571,25 @@ BOOST_AUTO_TEST_CASE( grid_dist_id_domain_test_use)
// Domain
Box<3,float> domain3({-0.3,-0.3,-0.3},{1.1,1.1,1.1});
periodicity<3> np({NON_PERIODIC,NON_PERIODIC,NON_PERIODIC});
periodicity<3> p({PERIODIC,PERIODIC,PERIODIC});
long int k = 128*128*128*create_vcluster().getProcessingUnits();
k = std::pow(k, 1/3.);
Test3D_domain(domain3,k);
Test3D_domain(domain3,k,np);
auto & v_cl = create_vcluster();
if (v_cl.getProcessingUnits() > 32)
return;
// We use a 128x128x128 and we move tha domain
for (size_t i = 0 ; i < 10 ; i++)
{
Box<3,float> exp({0.0,0.0,0.0},{1.3,1.3,1.3});
domain3.enlarge(exp);
Test3D_domain(domain3,128,p);
}
}
BOOST_AUTO_TEST_CASE( grid_dist_id_extended )
......
src/Grid/staggered_dist_grid.hpp
View file @ ac799346
......@@ -146,7 +146,7 @@ public:
auto g_dst_it = g_dst.getDomainIterator();
// Check that the 2 iterator has the same size
checkIterator<St,decltype(g_map_it),decltype(g_dst_it)>(g_map_it,g_dst_it);
checkIterator<Grid_dst,decltype(g_map_it),decltype(g_dst_it)>(g_map_it,g_dst_it);
while (g_map_it.isNext() == true)
{
......
src/Grid/staggered_dist_grid_copy.hpp
View file @ ac799346
......@@ -124,8 +124,8 @@ struct interp_ele
* this mean that the object passed must not be a temporal object
*
*/
inline interp_ele(grid_dist_key_dx<scheme::Sys_eqs_typ::dims> & key, S & grid_dst, Ev && x)
:key(key),grid_dst(grid_dst),x(x),interp_pos(interp_pos)
inline interp_ele(const grid_dist_key_dx<Grid_dst::dims> & key_dst, Grid_dst && grid_dst, const Grid_src & x, const grid_dist_key_dx<Grid_src