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
    struct point_and_gid
    {
        // global id of the particle
        size_t id;
 
        // Position of the neighborhood particle
        Point<3,float> xq;
 
        // Used to reorder the neighborhood particles by id
        bool operator<(const struct point_and_gid & pag) const
        {
            return (id < pag.id);
        }
    };
 

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
    typedef  aggregate<size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;

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;
    }
 

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;
    }
 

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>();
 

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
    struct point_and_gid
    {
        // global id of the particle
        size_t id;
 
        // Position of the neighborhood particle
        Point<3,float> xq;
 
        // Used to reorder the neighborhood particles by id
        bool operator<(const struct point_and_gid & pag) const
        {
            return (id < pag.id);
        }
    };
 
 
 
    // Particle properties list
    typedef  aggregate<size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
 
 
 
    // 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();
 
 
}