Commit 47964047 authored by Pietro Incardona's avatar Pietro Incardona

Semiworking periodic boundary conditions

parent 6e8741cf
openfpm_data @ 0572e97f
Subproject commit fb93d17019fdfb1f7c65c463143f5100d506d6d2
Subproject commit 0572e97f47cfe83ca1b3f66e05c47cda85946345
openfpm_devices @ 52cee143
Subproject commit ec5abaf0d363649dfa1d5831ac6784868ce418f2
Subproject commit 52cee143052e1b6ce4bdb0d6061d0e8095829643
openfpm_io @ 7988298e
Subproject commit c54eca6ee2e0bbfefda778df50b369ce9a90323f
Subproject commit 7988298e26138d0470851e9053e6640511793363
openfpm_vcluster @ e4f99633
Subproject commit d052b0d75a49835b24419cc4938485d6dc514393
Subproject commit e4f9963378be88d411b56d07afbeed05df994d9b
......@@ -489,6 +489,48 @@ public:
}
};
/*! \brief Apply boundary condition to the point
*
* \param p Point to apply the boundary condition
*
*/
void applyPointBC(float (& pt)[dim]) const
{
for (size_t i = 0 ; i < dim ; i++)
{
if (bc[i] == PERIODIC)
pt[i] = openfpm::math::periodic_l(pt[i],domain.getHigh(i),domain.getLow(i));
}
}
/*! \brief Apply boundary condition to the point
*
* \param p Point to apply the boundary condition
*
*/
void applyPointBC(Point<dim,T> & pt) const
{
for (size_t i = 0 ; i < dim ; i++)
{
if (bc[i] == PERIODIC)
pt.get(i) = openfpm::math::periodic_l(pt.get(i),domain.getHigh(i),domain.getLow(i));
}
}
/*! \brief Apply boundary condition to the point
*
* \param encapsulated object
*
*/
template<typename Mem> void applyPointBC(encapc<1,Point<dim,T>,Mem> && pt) const
{
for (size_t i = 0 ; i < dim ; i++)
{
if (bc[i] == PERIODIC)
pt.template get<0>()[i] = openfpm::math::periodic_l(pt.template get<0>()[i],domain.getHigh(i),domain.getLow(i));
}
}
/*! It calculate the internal ghost boxes
*
* Example: Processor 10 calculate
......@@ -789,6 +831,51 @@ p1[0]<-----+ +----> p2[0]
return fine_s.get(cd.getCell(p));
}
/*! \brief Given a point return in which processor the particle should go
*
* Boundary conditions are considered
*
* \return processorID
*
*/
template<typename Mem> size_t inline processorIDBC(encapc<1, Point<dim,T>, Mem> p)
{
Point<dim,T> pt = p;
applyPointBC(pt);
return fine_s.get(cd.getCell(pt));
}
/*! \brief Given a point return in which processor the particle should go
*
* Boundary conditions are considered
*
* \return processorID
*
*/
size_t inline processorIDBC(const Point<dim,T> &p) const
{
Point<dim,T> pt = p;
applyPointBC(pt);
return fine_s.get(cd.getCell(p));
}
/*! \brief Given a point return in which processor the particle should go
*
* Boundary consition are considered
*
* \return processorID
*
*/
size_t inline processorIDBC(const T (&p)[dim]) const
{
Point<dim,T> pt = p;
applyPointBC(pt);
return fine_s.get(cd.getCell(p));
}
/*! \brief Get the smallest subdivision of the domain on each direction
*
* \return a box p1 is set to zero
......@@ -879,7 +966,7 @@ p1[0]<-----+ +----> p2[0]
* \return The physical domain
*
*/
Domain<dim,T> & getDomain()
const Domain<dim,T> & getDomain()
{
return domain;
}
......
......@@ -191,7 +191,8 @@ struct N_box
openfpm::vector<comb<dim>> pos;
//! Default constructor
N_box() {};
N_box()
{};
//! Copy constructor
N_box(const N_box<dim,T> & b)
......@@ -213,8 +214,8 @@ struct N_box
N_box<dim,T> & operator=(const N_box<dim,T> & ele)
{
id = ele.id;
bx = ele.bx;
pos = ele.pos;
return * this;
}
......@@ -227,8 +228,8 @@ struct N_box
N_box<dim,T> & operator=(N_box<dim,T> && ele)
{
id = ele.id;
bx.swap(ele.bx);
pos = ele.pos;
return * this;
}
......@@ -243,6 +244,9 @@ struct N_box
if (id != ele.id)
return false;
if (pos != ele.pos)
return false;
return bx == ele.bx;
}
......
......@@ -156,13 +156,13 @@ protected:
switch (cmbs[j][k])
{
case 1:
shifts.get(cmbs[j].lin()).template get<0>()[0] = -domain.getHigh(k);
shifts.get(cmbs[j].lin()).template get<0>()[k] = -domain.getHigh(k);
break;
case 0:
shifts.get(cmbs[j].lin()).template get<0>()[0] = 0;
shifts.get(cmbs[j].lin()).template get<0>()[k] = 0;
break;
case -1:
shifts.get(cmbs[j].lin()).template get<0>()[0] = domain.getHigh(k);
shifts.get(cmbs[j].lin()).template get<0>()[k] = domain.getHigh(k);
break;
}
}
......@@ -395,7 +395,7 @@ public:
/*! It return the shift vector
*
* Consider a domain with some ghost, at the border of the domain the
* ghost must be threated in a special way depending on the periodicity
* ghost must be treated in a special way, depending on the periodicity
* of the boundary
*
\verbatim
......@@ -637,7 +637,7 @@ public:
* can produce more entry with the same processor, the UNIQUE option eliminate double entries
* (UNIQUE) is for particle data (MULTIPLE) is for grid data [default MULTIPLE]
*
* \param return the processor ids
* \param return the processor ids (not the rank, the id in the near processor list)
*
*/
template <typename id1, typename id2> inline const openfpm::vector<std::pair<size_t,size_t>> ghost_processorID_pair(Point<dim,T> & p, const int opt = MULTIPLE)
......
......@@ -282,6 +282,16 @@ public:
return *this;
}
/*! \brief Get the number of local sub-domains
*
* \return the number of local sub-domains
*
*/
inline size_t getNLocalSub()
{
return loc_ghost_box.size();
}
/*! \brief Get the number of external local ghost box for each sub-domain
*
* \param id sub-domain id
......@@ -337,6 +347,61 @@ public:
return loc_ghost_box.get(i).ibx.get(j).bx;
}
/*! \brief Get the j internal local ghost box boundary position for the i sub-domain of the local processor
*
* \note For the sub-domain i intersected with the sub-domain j enlarged, the associated
* external ghost box is located in getLocalIGhostBox(j,k) with
* getLocalIGhostSub(j,k) == i
*
* To get k use getLocalIGhostE
*
* \see getLocalIGhostE
*
* Some of the intersection boxes has special position, because they are at the boundary, this function
* return their position at the border
*
\verbatim
[1,1]
+---------+------------------------+---------+
| (1,-1) | | (1,1) |
| | | (1,0) --> 7 | | |
| v | | v |
| 6 | | 8 |
+--------------------------------------------+
| | | |
| | | |
| | | |
| (-1,0) | | (1,0) |
| | | | | |
| v | (0,0) --> 4 | v |
| 3 | | 5 |
| | | |
| | | |
| | | |
| | | |
| | | |
| | | |
+--------------------------------------------+
| (-1,-1) | | (-1,1) |
| | | (-1,0) --> 1 | | |
| v | | v |
| 0 | | 2 |
+---------+------------------------+---------+
\endverbatim
*
* \param i sub-domain
* \param j box
* \return the box
*
*/
inline const comb<dim> & getLocalIGhostPos(size_t i, size_t j) const
{
return loc_ghost_box.get(i).ibx.get(j).cmb;
}
/*! \brief Get the j external local ghost box for the local processor
*
* \param i sub-domain
......
......@@ -66,7 +66,6 @@ class nn_prcs
// cast the pointer
nn_prcs<dim,T> * cd = static_cast< nn_prcs<dim,T> *>(ptr);
// Resize the memory
cd->nn_processor_subdomains[i].bx.resize(msg_i / sizeof(::Box<dim,T>) );
// Return the receive pointer
......@@ -260,7 +259,7 @@ public:
{
for (size_t p = 0 ; p < getNNProcessors() ; p++)
{
auto list_p_box = getNearSubdomains(IDtoProc(p));
const openfpm::vector< ::Box<dim,T> > & list_p_box = getNearSubdomains(IDtoProc(p));
// Create the smallest box contained in all sub-domain
for (size_t b = 0 ; b < list_p_box.size() ; b++)
......@@ -328,6 +327,8 @@ public:
}
}
nn_processor_subdomains.reserve(nn_processors.size());
// Get the sub-domains of the near processors
v_cl.sendrecvMultipleMessagesNBX(nn_processors,boxes,nn_prcs<dim,T>::message_alloc, this ,NEED_ALL_SIZE);
......@@ -386,6 +387,7 @@ public:
std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " error this process rank is not adjacent to the local processor";
}
#endif
return key->second.bx;
}
......
......@@ -31,8 +31,6 @@
#define GET 1
#define PUT 2
#define INTERNAL 0
#define NO_POSITION 1
#define WITH_POSITION 2
......@@ -78,15 +76,354 @@ private:
//! Particle position vector, (It has 2 elements) the first has real particles assigned to a processor
//! the second element contain unassigned particles
Vcluster_object_array<openfpm::vector<Point<dim,St>>> v_pos;
openfpm::vector<Point<dim,St>> v_pos;
//! Particle properties vector, (It has 2 elements) the first has real particles assigned to a processor
//! the second element contain unassigned particles
Vcluster_object_array<openfpm::vector<prop>> v_prp;
openfpm::vector<prop> v_prp;
//! Virtual cluster
Vcluster & v_cl;
// definition of the send vector for position
typedef openfpm::vector<Point<dim,St>,ExtPreAlloc<Memory>> send_pos_vector;
//////////////////////////////
// COMMUNICATION variables
//////////////////////////////
//! It map the processor id with the communication request into map procedure
openfpm::vector<size_t> p_map_req;
//! For each near processor, outgoing particle id and shift vector
openfpm::vector<openfpm::vector<std::pair<size_t,size_t>>> opart;
//! For each adjacent processor the size of the ghost sending buffer
openfpm::vector<size_t> ghost_prc_sz;
//! Sending buffer for the ghost particles properties
Memory g_prp_mem;
//! Sending buffer for the ghost particles position
Memory g_pos_mem;
//! For each adjacent processor it store the size of the receiving message in byte
openfpm::vector<size_t> recv_sz;
//! For each adjacent processot it store the received message
openfpm::vector<HeapMemory> recv_mem_gg;
//! Receive buffer for global communication
HeapMemory hp_recv;
//! For each message contain the processor from which processor come from
openfpm::vector<size_t> v_proc;
//! Total size of the received buffer
size_t recv_cnt;
/*! \brief Label the particles
*
* It count the number of particle to send to each processors and save its ids
*
*/
void labelParticles()
{
// Buffer that contain the number of elements to send for each processor
ghost_prc_sz.clear();
ghost_prc_sz.resize(dec.getNNProcessors());
// Buffer that contain for each processor the id of the particle to send
opart.clear();
opart.resize(dec.getNNProcessors());
// Iterate over all particles
auto it = v_pos.getIterator();
while (it.isNext())
{
auto key = it.get();
// Given a particle, it return which processor require it (first id) and shift id, second id
// For an explanation about shifts vectors please consult getShiftVector in ie_ghost
const openfpm::vector<std::pair<size_t,size_t>> & vp_id = dec.template ghost_processorID_pair<typename Decomposition::lc_processor_id, typename Decomposition::shift_id>(v_pos.get(key),UNIQUE);
for (size_t i = 0 ; i < vp_id.size() ; i++)
{
// processor id
size_t p_id = vp_id.get(i).first;
// add particle to communicate
ghost_prc_sz.get(p_id)++;
opart.get(p_id).add(std::pair<size_t,size_t>(key,vp_id.get(i).second));
}
++it;
}
}
/*! \brief Add local particles based on the boundary conditions
*
* In order to understand what this function use the following
*
\verbatim
[1,1]
+---------+------------------------+---------+
| (1,-1) | | (1,1) |
| | | (1,0) --> 7 | | |
| v | | v |
| 6 | | 8 |
+--------------------------------------------+
| | | |
| | | |
| | | |
| (-1,0) | | (1,0) |
| | | | | |
| v | (0,0) --> 4 | v |
| 3 | | 5 |
| | | |
B | | | A |
* | | | * |
| | | |
| | | |
| | | |
+--------------------------------------------+
| (-1,-1) | | (-1,1) |
| | | (-1,0) --> 1 | | |
| v | | v |
| 0 | | 2 |
+---------+------------------------+---------+
\endverbatim
*
* The box is the domain, while all boxes at the border (so not (0,0) ) are the
* ghost part at the border of the domain. If a particle A is in the position in figure
* a particle B must be created. This function duplicate the particle A, if A and B are
* local
*
*/
void add_loc_particles_bc()
{
// get the shift vectors
const openfpm::vector<Point<dim,St>> & shifts = dec.getShiftVectors();
// Add local particles coming from periodic boundary, the only boxes that count are the one
// at the border, filter them
openfpm::vector_std<Box<dim,St>> box_f;
openfpm::vector_std<comb<dim>> box_cmb;
for (size_t i = 0 ; i < dec.getNLocalSub() ; i++)
{
size_t Nl = dec.getLocalNIGhost(i);
for (size_t j = 0 ; j < Nl ; j++)
{
// If they are not in the border the combination is all zero
if (dec.getLocalIGhostPos(i,j).n_zero() == dim)
continue;
box_f.add(dec.getLocalIGhostBox(i,j));
box_cmb.add(dec.getLocalIGhostPos(i,j));
}
}
if (box_f.size() == 0)
return;
else
{
// Label the internal (assigned) particles
auto it = v_pos.getIterator();
while (it.isNext())
{
auto key = it.get();
// 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)
{
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);
}
}
++it;
}
}
}
/*! \brief This function fill the send buffer for the particle position after the particles has been label with labelParticles
*
* \param g_pos_send Send buffer to fill
* \param prAlloc_pos Memory object for the send buffer
*
*/
void fill_send_pos_buf(openfpm::vector<send_pos_vector> & g_pos_send, ExtPreAlloc<Memory> * prAlloc_pos)
{
// get the shift vectors
const openfpm::vector<Point<dim,St>> & shifts = dec.getShiftVectors();
// create a number of send buffers equal to the near processors
g_pos_send.resize(ghost_prc_sz.size());
for (size_t i = 0 ; i < g_pos_send.size() ; i++)
{
// set the preallocated memory to ensure contiguity
g_pos_send.get(i).setMemory(*prAlloc_pos);
// resize the sending vector (No allocation is produced)
g_pos_send.get(i).resize(ghost_prc_sz.get(i));
}
// Fill the send buffer
for ( size_t i = 0 ; i < opart.size() ; i++ )
{
for (size_t j = 0 ; j < opart.get(i).size() ; j++)
{
Point<dim,St> s = v_pos.get(opart.get(i).get(j).first);
s -= shifts.get(opart.get(i).get(j).second);
g_pos_send.get(i).set(j,s);
}
}
}
/*! \brief This function fill the send buffer for properties after the particles has been label with labelParticles
*
* \tparam send_vector type used to send data
* \tparam prp_object object containing only the properties to send
* \tparam prp set of properties to send
*
* \param g_send_prp Send buffer to fill
* \param prAlloc_prp Memory object for the send buffer
*
*/
template<typename send_vector,typename prp_object, int... prp> void fill_send_prp_buf(openfpm::vector<send_vector> & g_send_prp, ExtPreAlloc<Memory> * prAlloc_prp)
{
// create a number of send buffers equal to the near processors
g_send_prp.resize(ghost_prc_sz.size());
for (size_t i = 0 ; i < g_send_prp.size() ; i++)
{
// set the preallocated memory to ensure contiguity
g_send_prp.get(i).setMemory(*prAlloc_prp);
// resize the sending vector (No allocation is produced)
g_send_prp.get(i).resize(ghost_prc_sz.get(i));
}
// Fill the send buffer
for ( size_t i = 0 ; i < opart.size() ; i++ )
{
for (size_t j = 0 ; j < opart.get(i).size() ; j++)
{
// source object type
typedef encapc<1,prop,typename openfpm::vector<prop>::memory_conf> encap_src;
// destination object type
typedef encapc<1,prp_object,typename openfpm::vector<prp_object>::memory_conf> encap_dst;
// Copy only the selected properties
object_si_d<encap_src,encap_dst,OBJ_ENCAP,prp...>(v_prp.get(opart.get(i).get(j).first),g_send_prp.get(i).get(j));
}
}
}
/*! \brief This function process the receiced data for the properties and populate the ghost
*
* \tparam send_vector type used to send data
* \tparam prp_object object containing only the properties to send
* \tparam prp set of properties to send
*
*/
template<typename send_vector,typename prp_object, int... prp> void process_received_prp()
{
// Mark the ghost part
g_m = v_prp.size();
// Process the received data (recv_mem_gg != 0 if you have data)
for (size_t i = 0 ; i < dec.getNNProcessors() && recv_mem_gg.size() != 0 ; i++)
{
// calculate the number of received elements
size_t n_ele = recv_sz.get(i) / sizeof(prp_object);
// add the received particles to the vector
PtrMemory * ptr1 = new PtrMemory(recv_mem_gg.get(i).getPointer(),recv_sz.get(i));
// create vector representation to a piece of memory already allocated
openfpm::vector<prp_object,PtrMemory,openfpm::grow_policy_identity> v2;
v2.setMemory(*ptr1);
// resize with the number of elements
v2.resize(n_ele);
// Add the ghost particle
v_prp.template add_prp<prp_object,PtrMemory,openfpm::grow_policy_identity,prp...>(v2);
}
}
/*! \brief This function process the received data for the properties and populate the ghost
*
*/
void process_received_pos()
{
// Process the received data (recv_mem_gg != 0 if you have data)
for (size_t i = 0 ; i < dec.getNNProcessors() && recv_mem_gg.size() != 0 ; i++)
{
// calculate the number of received elements
size_t n_ele = recv_sz.get(i) / sizeof(Point<dim,St>);
// add the received particles to the vector
PtrMemory * ptr1 = new PtrMemory(recv_mem_gg.get(i).getPointer(),recv_sz.get(i));
// create vector representation to a piece of memory already allocated
openfpm::vector<Point<dim,St>,PtrMemory,openfpm::grow_policy_identity> v2;
v2.setMemory(*ptr1);
// resize with the number of elements
v2.resize(n_ele);