OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
Vector 6 complex usage for validation and debugging

Complex usage for validation and debugging

In this example we show how the flexibility of the library can be used to perform complex tasks for validation and debugging. We will use a lot of feature that has been explained in the previous examples

In a previous example we show in case of symmetric interaction between particles how to symmetric cell list could be used to speed up the calculation.

See also
not present

In this example we validate that both computation match, in particular because the computation depend from the neighborhood, we will check and validate that the neighborhood of each particle is equivalent with symmetric cell list and normal cell-list.

Initialization

The initialization is classical we define cutoff radius, domain, boundary conditions, ghost part and some useful constants. We also define a struct that will contain the debug information. This structure contain an id that is the global id of the particle and a point that is the the position of the neighborhood.

See also
Initialization
openfpm_init(&argc,&argv);
constexpr int gid = 0;
constexpr int nn_norm = 1;
constexpr int nn_sym = 2;
// used to define the domain
float L = 1000.0;
// Domain
Box<3,float> box({-L,-L,-L},{L,L,L});
// Boundary conditions
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
// cut-off radius
float r_cut = 100.0;
// ghost
Ghost<3,float> ghost(r_cut);
// Point and global id
{
// global id of the particle
size_t id;
// Position of the neighborhood particle
// Used to reorder the neighborhood particles by id
bool operator<(const struct point_and_gid & pag) const
{
return (id < pag.id);
}
};
This class represent an N-dimensional box.
Definition Box.hpp:61
This class implement the point shape in an N-dimensional space.
Definition Point.hpp:28
Definition ids.hpp:149

Particle

For this example we will use a particle with 3 properties

  • The first property is a global index for the particle unique across processors
  • The second is a vector that contain for each neighborhood the global id of the particle and its position
  • The last one is again a vector equivalent to the previous but is produced with the symmetric cell-list
// Particle properties list
Implementation of 1-D std::vector like structure.
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...

Distributed vector

Here we create a distributed vector of 4096 particles. We initialize the particles randomly and we assign a global index to the first property of the particle

The function accum return the following number

\( \sum_{i < proc} size\_local(i) \)

where \( proc \) is the processor where we are computing the formula and \( size\_local(i) \) is the number of particles the processor \( i \) has. This number is used to produce a global id of the particles

// Distributed vector
vector_dist<3,float, part_prop > vd(4096,box,bc,ghost);
// used to calculate a global index for each particle
size_t start = vd.accum();
auto it = vd.getDomainIterator();
while (it.isNext())
{
auto key = it.get();
vd.getPos(key)[0] = 2.0*L*((float)rand()/RAND_MAX) - L;
vd.getPos(key)[1] = 2.0*L*((float)rand()/RAND_MAX) - L;
vd.getPos(key)[2] = 2.0*L*((float)rand()/RAND_MAX) - L;
vd.getProp<gid>(key) = key.getKey() + start;
++it;
}
Distributed vector.

Redistribute

Redistribute the particles and synchronize the ghosts. In this case we are only interested in synchronizing the global id property.

vd.map();
// sync the ghost
vd.ghost_get<gid>();

Calculate the neighborhood of each particle

Here we calculate the neighborhood of each particles using a simple cell list. In case the particle q has a distance from the particle p smaller than r_cut. We add it in the first list of neighborhood

auto NN = vd.getCellList(r_cut);
auto p_it = vd.getDomainIterator();
while (p_it.isNext())
{
auto p = p_it.get();
Point<3,float> xp = vd.getPos(p);
auto Np = NN.getNNIterator(NN.getCell(vd.getPos(p)));
while (Np.isNext())
{
auto q = Np.get();
// skip self interaction
if (p.getKey() == q)
{
++Np;
continue;
}
Point<3,float> xq = vd.getPos(q);
Point<3,float> f = (xp - xq);
float distance = f.norm();
// if the distance smalle than the cut-off radius add it to the neighborhood list
if (distance < r_cut )
{
vd.getProp<nn_norm>(p).add();
vd.getProp<nn_norm>(p).last().xq = xq;
vd.getProp<nn_norm>(p).last().id = vd.getProp<0>(q);
}
// Next neighborhood
++Np;
}
// Next particle
++p_it;
}
__device__ __host__ const T & get(unsigned int i) const
Get coordinate.
Definition Point.hpp:172

Calculate the neighborhood of each particle with symmetric cell-list

Here we calculate the neighborhood of each particle using instead a symmetric cell list. In case of symmetric cell-list if we find that a particle p is neighborhood of particle q we have to add p to the neighborhood of q and q to the neighborhood of p. Because q can be a ghost particle when we add the neighborhood p to q we have to transmit such information to the real owner of the particle. This can be done with the function ghost_put. In this case we use the operation merge_ that add the already calculated neighborhood with the transmitted one. More in general it merge the information together.

auto NN2 = vd.getCellListSym(r_cut);
auto p_it2 = vd.getDomainIterator();
while (p_it2.isNext())
{
auto p = p_it2.get();
Point<3,float> xp = vd.getPos(p);
auto Np = NN2.template getNNIteratorSym<NO_CHECK>(NN2.getCell(vd.getPos(p)),p.getKey(),vd.getPosVector());
while (Np.isNext())
{
auto q = Np.get();
if (p.getKey() == q)
{
++Np;
continue;
}
// 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<nn_sym>(p).add();
vd.getProp<nn_sym>(q).add();
vd.getProp<nn_sym>(p).last().xq = xq;
vd.getProp<nn_sym>(q).last().xq = xp;
vd.getProp<nn_sym>(p).last().id = vd.getProp<0>(q);
vd.getProp<nn_sym>(q).last().id = vd.getProp<0>(p);
}
++Np;
}
++p_it2;
}
vd.ghost_put<merge_,2>();
This structure define the operation add to use with copy general.

Cheking for validation

Here we check that the two calculated neighborhood match. In particular, because the order of the particles does not match, we have to first reorder by global-id, and than check the list. We cannot instead easily compare the position. The reason can be seen in this figure

     +----------------+
     |                |
     |           1    |
     |           *  2 |
     |           3  * |
   +<----- ghost *  +<-------- real
     |* <--- real     |* <------- ghost
     |4               |4
     |                |
     |                |
     +----------------+

In the case we are calculating the neighborhood of the particle \( + \) in case of normal Cell-list. The real particle 1,2,3 are added, while the particle 4 is added with the ghost particle coordinates. In the case of symmetric cell-list the real 1,2,3 are added to \( + \). The particle 4 instead is added by the ghost put. In particular 4 real add the particle to the \( + \) ghost particle and than ghost put merge the information. This mean that the symmetric add the particle 4 with the real coordinates of the particle 4

auto p_it3 = vd.getDomainIterator();
bool ret = true;
while (p_it3.isNext())
{
auto p = p_it3.get();
vd.getProp<nn_norm>(p).sort();
vd.getProp<nn_sym>(p).sort();
ret &= vd.getProp<nn_norm>(p).size() == vd.getProp<nn_sym>(p).size();
for (size_t i = 0 ; i < vd.getProp<1>(p).size() ; i++)
{
ret &= vd.getProp<nn_norm>(p).get(i).id == vd.getProp<nn_sym>(p).get(i).id;
if (box.isInside(vd.getProp<nn_norm>(p).get(i).xq) == true)
{
ret &= vd.getProp<nn_norm>(p).get(i).xq == vd.getProp<nn_sym>(p).get(i).xq;
}
}
++p_it3;
}
if (ret != true)
{
std::cout << "ERROR" << std::endl;
exit(1);
}

Finalize

At the very end of the program we have always to de-initialize the library

openfpm_finalize();

Full code

#include "Vector/vector_dist.hpp"
#include "data_type/aggregate.hpp"
#include "Plot/GoogleChart.hpp"
#include "Plot/util.hpp"
#include "timer.hpp"
int main(int argc, char* argv[])
{
openfpm_init(&argc,&argv);
constexpr int gid = 0;
constexpr int nn_norm = 1;
constexpr int nn_sym = 2;
// used to define the domain
float L = 1000.0;
// Domain
Box<3,float> box({-L,-L,-L},{L,L,L});
// Boundary conditions
size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
// cut-off radius
float r_cut = 100.0;
// ghost
Ghost<3,float> ghost(r_cut);
// Point and global id
{
// global id of the particle
size_t id;
// Position of the neighborhood particle
// Used to reorder the neighborhood particles by id
bool operator<(const struct point_and_gid & pag) const
{
return (id < pag.id);
}
};
// Particle properties list
// Distributed vector
vector_dist<3,float, part_prop > vd(4096,box,bc,ghost);
// used to calculate a global index for each particle
size_t start = vd.accum();
auto it = vd.getDomainIterator();
while (it.isNext())
{
auto key = it.get();
vd.getPos(key)[0] = 2.0*L*((float)rand()/RAND_MAX) - L;
vd.getPos(key)[1] = 2.0*L*((float)rand()/RAND_MAX) - L;
vd.getPos(key)[2] = 2.0*L*((float)rand()/RAND_MAX) - L;
vd.getProp<gid>(key) = key.getKey() + start;
++it;
}
vd.map();
// sync the ghost
vd.ghost_get<gid>();
auto NN = vd.getCellList(r_cut);
auto p_it = vd.getDomainIterator();
while (p_it.isNext())
{
auto p = p_it.get();
Point<3,float> xp = vd.getPos(p);
auto Np = NN.getNNIterator(NN.getCell(vd.getPos(p)));
while (Np.isNext())
{
auto q = Np.get();
// skip self interaction
if (p.getKey() == q)
{
++Np;
continue;
}
Point<3,float> xq = vd.getPos(q);
Point<3,float> f = (xp - xq);
float distance = f.norm();
// if the distance smalle than the cut-off radius add it to the neighborhood list
if (distance < r_cut )
{
vd.getProp<nn_norm>(p).add();
vd.getProp<nn_norm>(p).last().xq = xq;
vd.getProp<nn_norm>(p).last().id = vd.getProp<0>(q);
}
// Next neighborhood
++Np;
}
// Next particle
++p_it;
}
auto NN2 = vd.getCellListSym(r_cut);
auto p_it2 = vd.getDomainIterator();
while (p_it2.isNext())
{
auto p = p_it2.get();
Point<3,float> xp = vd.getPos(p);
auto Np = NN2.template getNNIteratorSym<NO_CHECK>(NN2.getCell(vd.getPos(p)),p.getKey(),vd.getPosVector());
while (Np.isNext())
{
auto q = Np.get();
if (p.getKey() == q)
{
++Np;
continue;
}
// 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<nn_sym>(p).add();
vd.getProp<nn_sym>(q).add();
vd.getProp<nn_sym>(p).last().xq = xq;
vd.getProp<nn_sym>(q).last().xq = xp;
vd.getProp<nn_sym>(p).last().id = vd.getProp<0>(q);
vd.getProp<nn_sym>(q).last().id = vd.getProp<0>(p);
}
++Np;
}
++p_it2;
}
vd.ghost_put<merge_,2>();
auto p_it3 = vd.getDomainIterator();
bool ret = true;
while (p_it3.isNext())
{
auto p = p_it3.get();
vd.getProp<nn_norm>(p).sort();
vd.getProp<nn_sym>(p).sort();
ret &= vd.getProp<nn_norm>(p).size() == vd.getProp<nn_sym>(p).size();
for (size_t i = 0 ; i < vd.getProp<1>(p).size() ; i++)
{
ret &= vd.getProp<nn_norm>(p).get(i).id == vd.getProp<nn_sym>(p).get(i).id;
if (box.isInside(vd.getProp<nn_norm>(p).get(i).xq) == true)
{
ret &= vd.getProp<nn_norm>(p).get(i).xq == vd.getProp<nn_sym>(p).get(i).xq;
}
}
++p_it3;
}
if (ret != true)
{
std::cout << "ERROR" << std::endl;
exit(1);
}
openfpm_finalize();
}