vector_dist_gpu_MP_tests.cu

#define BOOST_TEST_DYN_LINK
#include <boost/test/unit_test.hpp>
 
#include "Vector/vector_dist_multiphase_functions.hpp"
#include "VCluster/VCluster.hpp"
#include "Vector/vector_dist.hpp"
 
BOOST_AUTO_TEST_SUITE( vector_dist_multiphase_gpu_test )
 
BOOST_AUTO_TEST_CASE( vector_dist_multiphase_gpu_cell_list_test )
{
    if (create_vcluster().getProcessingUnits() > 24)
        return;
 
    size_t sz[3] = {60,60,40};
 
    // The domain
    Box<3,float> box({-1000.0,-1000.0,-1000.0},{2000.0,2000.0,1000.0});
 
    // Boundary conditions
    size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
 
    float r_cut = 51.0;
 
    // ghost, big enough to contain the interaction radius
    Ghost<3,float> ghost(r_cut);
 
    openfpm::vector< vector_dist_gpu<3,float, aggregate<double,double>> > phases;
 
    // first phase
    phases.add( vector_dist_gpu<3,float, aggregate<double,double>>(0,box,bc,ghost) );
 
    // The other 3 phases
    phases.add( vector_dist_gpu<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),0) );
    phases.add( vector_dist_gpu<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),0) );
    phases.add( vector_dist_gpu<3,float, aggregate<double,double>>(phases.get(0).getDecomposition(),0) );
 
    // Fill the phases with particles
 
    auto g_it = phases.get(0).getGridIterator(sz);
 
    while (g_it.isNext())
    {
        auto key = g_it.get();
 
        // Add a particle to all the phases
        phases.get(0).add();
        phases.get(1).add();
        phases.get(2).add();
        phases.get(3).add();
 
        phases.get(0).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
        phases.get(0).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
        phases.get(0).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
 
        phases.get(1).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
        phases.get(1).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
        phases.get(1).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
 
        phases.get(2).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
        phases.get(2).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
        phases.get(2).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
 
        phases.get(3).getLastPos()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
        phases.get(3).getLastPos()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
        phases.get(3).getLastPos()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
 
        ++g_it;
    }
 
    // Sync all phases
    for (size_t i = 0 ; i < 4 ; i++)
    {
        phases.get(i).map();
    }
 
    // randomize a little the particles
 
    for (size_t p = 0 ; p < phases.size() ; p++)
    {
        openfpm::vector_gpu<Point<3,float>> vt;
 
        for (size_t j = 0 ; j < phases.get(p).size_local() ; j++)
        {
            vt.add(phases.get(p).getPos((j + p*133) % phases.get(p).size_local()));
        }
        phases.get(p).getPosVector().swap(vt);
    }
 
    // Sync all phases
    for (size_t i = 0 ; i < 4 ; i++)
    {
        phases.get(i).ghost_get<>();
    }
 
    // Get the cell list of the phase 0 and 1
    auto CL_phase0 = phases.get(0).getCellList(r_cut);
    auto CL_phase1 = phases.get(1).getCellList(r_cut);
 
    // This function create a Verlet-list between phases 0 and 1
    auto NN_ver01 = createVerlet(phases.get(0),phases.get(1),CL_phase1,r_cut);
 
    // Check NNver0_1
 
    bool ret = true;
    auto it = phases.get(0).getDomainIterator();
 
    while (it.isNext())
    {
        auto p = it.get();
        auto Np = NN_ver01.getNNIterator<NO_CHECK>(p.getKey());
 
        size_t nn_count = 0;
 
        // For each neighborhood of the particle p
        while (Np.isNext())
        {
            // Count the number of particles
            nn_count++;
 
            ++Np;
        }
 
        ret &= nn_count == 7ul;
 
        ++it;
    }
 
    BOOST_REQUIRE_EQUAL(ret,true);
 
    // Sync all phases
    for (size_t i = 0 ; i < 4 ; i++)
    {
        phases.get(i).map();
        phases.get(i).ghost_get<>();
    }
 
    // NN_ver0_all
 
    // This function create an "Empty" Multiphase Cell List
    auto CL_all = createCellListM<2>(phases,r_cut);
 
    // This create a Verlet-list between phase 0 and all the other phases
    auto NNver0_all = createVerletM<2>(0,phases.get(0),phases,CL_all,r_cut);
 
    it = phases.get(0).getDomainIterator();
 
    while (it.isNext())
    {
        auto p = it.get();
        auto Np = NNver0_all.getNNIterator<NO_CHECK>(p.getKey());
 
        size_t nn_cout[4] = {0,0,0,0};
 
        // For each neighborhood of the particle p
        while (Np.isNext())
        {
            // Get from which phase it come from
            auto ph_q = Np.getV();
 
            nn_cout[ph_q]++;
 
            ++Np;
        }
 
        ret &= nn_cout[0] == 7;
        ret &= nn_cout[1] == 7;
        ret &= nn_cout[2] == 7;
        ret &= nn_cout[3] == 7;
 
        ++it;
    }
 
    BOOST_REQUIRE_EQUAL(ret,true);
}
 
 
BOOST_AUTO_TEST_CASE( vector_dist_multiphase_cell_list_sym_test )
{
    if (create_vcluster().getProcessingUnits() > 24)
        return;
 
    size_t sz[3] = {60,60,40};
 
    // The domain
    Box<3,float> box({-1000.0,-1000.0,-1000.0},{2000.0,2000.0,1000.0});
 
    // Boundary conditions
    size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
 
    float r_cut = 51.0;
 
    // ghost, big enough to contain the interaction radius
    Ghost<3,float> ghost(r_cut);
 
    openfpm::vector< vector_dist_gpu<3,float, aggregate<size_t>> > phases;
 
    // first phase
    phases.add( vector_dist_gpu<3,float, aggregate<size_t>>(0,box,bc,ghost) );
 
    // The other 3 phases
    phases.add( vector_dist_gpu<3,float, aggregate<size_t>>(phases.get(0).getDecomposition(),0) );
    phases.add( vector_dist_gpu<3,float, aggregate<size_t>>(phases.get(0).getDecomposition(),0) );
    phases.add( vector_dist_gpu<3,float, aggregate<size_t>>(phases.get(0).getDecomposition(),0) );
 
    // Fill the phases with particles
 
    auto g_it = phases.get(0).getGridIterator(sz);
 
    while (g_it.isNext())
    {
        auto key = g_it.get();
 
        // Add a particle to all the phases
        phases.get(0).add();
        phases.get(1).add();
        phases.get(2).add();
        phases.get(3).add();
 
        phases.get(0).getLastPosWrite()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
        phases.get(0).getLastPosWrite()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
        phases.get(0).getLastPosWrite()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
        phases.get(0).getLastPropWrite<0>() = 0;
 
        phases.get(1).getLastPosWrite()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
        phases.get(1).getLastPosWrite()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
        phases.get(1).getLastPosWrite()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
        phases.get(1).getLastPropWrite<0>() = 0;
 
        phases.get(2).getLastPosWrite()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
        phases.get(2).getLastPosWrite()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
        phases.get(2).getLastPosWrite()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
        phases.get(2).getLastPropWrite<0>() = 0;
 
        phases.get(3).getLastPosWrite()[0] = key.get(0) * g_it.getSpacing(0) + box.getLow(0);
        phases.get(3).getLastPosWrite()[1] = key.get(1) * g_it.getSpacing(1) + box.getLow(1);
        phases.get(3).getLastPosWrite()[2] = key.get(2) * g_it.getSpacing(2) + box.getLow(2);
        phases.get(3).getLastPropWrite<0>() = 0;
 
        ++g_it;
    }
 
    // Sync all phases
    for (size_t i = 0 ; i < 4 ; i++)
    {
        phases.get(i).map();
    }
 
    // randomize a little the particles
 
    for (size_t p = 0 ; p < phases.size() ; p++)
    {
        openfpm::vector_gpu<Point<3,float>> vt;
 
        for (size_t j = 0 ; j < phases.get(p).size_local() ; j++)
        {
            vt.add(phases.get(p).getPos((j + p*133) % phases.get(p).size_local()));
        }
        phases.get(p).getPosVector().swap(vt);
    }
 
    // Sync all phases
    for (size_t i = 0 ; i < 4 ; i++)
    {
        phases.get(i).ghost_get<0>();
    }
 
    // Get the cell list of the phase 0 and 1
    auto CL_phase0 = phases.get(0).getCellListSym(r_cut);
    auto CL_phase1 = phases.get(1).getCellListSym(r_cut);
 
    // This function create a Verlet-list between phases 0 and 1
    auto NN_ver01 = createVerletSym(phases.get(0),phases.get(1),CL_phase1,r_cut);
 
    // Check NNver0_1
 
    bool ret = true;
    auto it = phases.get(0).getDomainIterator();
 
    while (it.isNext())
    {
        auto p = it.get();
        auto Np = NN_ver01.getNNIterator<NO_CHECK>(p.getKey());
 
        // For each neighborhood of the particle p
        while (Np.isNext())
        {
            // Neighborhood particle q
            auto q = Np.get();
 
            phases.get(0).getPropWrite<0>(p)++;
            phases.get(1).getPropWrite<0>(q)++;
 
            ++Np;
        }
 
        ++it;
    }
 
    phases.get(0).ghost_put<add_,0>();
    phases.get(1).ghost_put<add_,0>();
 
#ifdef SE_CLASS3
 
    phases.get(1).getDomainIterator();
 
#endif
 
    it = phases.get(0).getDomainIterator();
    while (it.isNext())
    {
        auto p = it.get();
 
        ret &= phases.get(0).getPropRead<0>(p) == 7;
        ret &= phases.get(1).getPropRead<0>(p) == 7;
 
        ++it;
    }
 
    BOOST_REQUIRE_EQUAL(ret,true);
 
    // Sync all phases
    for (size_t i = 0 ; i < 4 ; i++)
    {
        phases.get(i).map();
        phases.get(i).ghost_get<>();
    }
 
    // Reset counter on all phases
 
    for (size_t i = 0 ; i < phases.size() ; i++)
    {
        it = phases.get(i).getDomainAndGhostIterator();
        while (it.isNext())
        {
            auto p = it.get();
 
            phases.get(i).getPropWrite<0>(p) = 0;
 
            ++it;
        }
    }
 
    // NN_ver0_all
 
    // This function create an "Empty" Multiphase Cell List
    auto CL_all = createCellListSymM<2>(phases,r_cut);
 
    typedef decltype(createVerletSymM<2>(0,phases.get(0),phases,CL_all,r_cut)) verlet_type;
 
    verlet_type NNver_all[4];
 
    // This create a Verlet-list between each phase to all the other phases
    NNver_all[0] = createVerletSymM<2>(0,phases.get(0),phases,CL_all,r_cut);
    NNver_all[1] = createVerletSymM<2>(1,phases.get(1),phases,CL_all,r_cut);
    NNver_all[2] = createVerletSymM<2>(2,phases.get(2),phases,CL_all,r_cut);
    NNver_all[3] = createVerletSymM<2>(3,phases.get(3),phases,CL_all,r_cut);
 
    // all phases to all phases
 
    for (size_t i = 0 ; i < phases.size() ; i++)
    {
        it = phases.get(i).getDomainIterator();
 
        while (it.isNext())
        {
            auto p = it.get();
            auto Np = NNver_all[i].getNNIterator<NO_CHECK>(p.getKey());
 
            // For each neighborhood of the particle p
            while (Np.isNext())
            {
                // Get the particle q near to p
                auto q = Np.getP();
 
                // Get from which phase it come from
                auto ph_q = Np.getV();
 
                phases.get(i).getPropWrite<0>(p)++;
                phases.get(ph_q).getPropWrite<0>(q)++;
 
                ++Np;
            }
 
            ++it;
        }
    }
 
    phases.get(0).ghost_put<add_,0>();
    phases.get(1).ghost_put<add_,0>();
    phases.get(2).ghost_put<add_,0>();
    phases.get(3).ghost_put<add_,0>();
 
#ifdef SE_CLASS3
 
    it = phases.get(1).getDomainIterator();
    it = phases.get(2).getDomainIterator();
    it = phases.get(3).getDomainIterator();
 
#endif
 
    it = phases.get(0).getDomainIterator();
    while (it.isNext())
    {
        auto p = it.get();
 
        ret &= phases.get(0).getPropRead<0>(p) == 32;
        ret &= phases.get(1).getPropRead<0>(p) == 32;
        ret &= phases.get(2).getPropRead<0>(p) == 32;
        ret &= phases.get(3).getPropRead<0>(p) == 32;
 
        if (ret == false)
        {
            std::cout << phases.get(0).getPropRead<0>(p) << std::endl;
            std::cout << phases.get(1).getPropRead<0>(p) << std::endl;
            std::cout << phases.get(2).getPropRead<0>(p) << std::endl;
            std::cout << phases.get(3).getPropRead<0>(p) << std::endl;
        }
 
        ++it;
    }
 
    BOOST_REQUIRE_EQUAL(ret,true);
}
 
template<typename mp_vector_type, typename output_type>
__global__ void vdmkt(mp_vector_type mp_v, output_type ot)
{
    int k = 0;
    for (int i = 0 ; i < mp_v.size() ; i++)
    {
        for (int j = 0 ; j < mp_v.template get<0>(i).size() ; j++)
        {
            ot.template get<0>(k) = mp_v.template get<0>(i). template getProp<0>(j);
            k++;
        }
    }
}
 
template<typename mp_vector_type, typename output_type>
__global__ void vdmkt_simple(mp_vector_type mp_v, output_type ot)
{
    int k = 0;
 
    for (int i = 0 ; i < mp_v.size() ; i++)
    {
        for (int j = 0 ; j < mp_v.get(i).size_local() ; j++)
        {
            ot.template get<0>(k) = mp_v.get(i).template getProp<0>(j);
            k++;
        }
    }
}
 
template<typename mp_vector_type, typename output_type, typename cl_type>
__global__ void vdmkt_simple_cl(mp_vector_type mp_v, output_type ot, cl_type cl, output_type ot2)
{
    int k = 0;
 
    for (int i = 0 ; i < mp_v.size() ; i++)
    {
        for (int j = 0 ; j < mp_v.get(i).size_local() ; j++)
        {
            ot.template get<0>(k) = mp_v.get(i).template getProp<0>(j);
            k++;
        }
    }
 
    k = 0;
    for (int i = 0 ; i < cl.size() ; i++)
    {
        for (int j = 0 ; j < cl.get(i).getNCells() ; j++)
        {
            for (int k = 0 ; k < cl.get(i).getNelements(j) ; k++)
            {
                auto s = cl.get(i).get(j,k);
 
                ot2.template get<0>(k) = s;
 
                k++;
            }
        }
    }
 
    //auto & cl_ph = cl.template get(0).getNNIterator();
}
 
BOOST_AUTO_TEST_CASE( vector_dist_multiphase_kernel_test )
{
    if (create_vcluster().getProcessingUnits() > 24)
    {return;}
 
    // The domain
    Box<3,float> box({-1000.0,-1000.0,-1000.0},{2000.0,2000.0,1000.0});
 
    // Boundary conditions
    size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
 
    float r_cut = 51.0;
 
    // ghost, big enough to contain the interaction radius
    Ghost<3,float> ghost(r_cut);
 
    openfpm::vector_gpu<aggregate<vector_dist_gpu<3,float, aggregate<float>>> > phases;
 
    // first phase
    phases.add_no_device();
    phases.template get<0>(0) = vector_dist_gpu<3,float, aggregate<float>>(0,box,bc,ghost);
 
    // The other 3 phases
    phases.add_no_device();
    phases.template get<0>(1) = vector_dist_gpu<3,float, aggregate<float>>(phases.template get<0>(0).getDecomposition(),0);
    phases.add_no_device();
    phases.template get<0>(2) = vector_dist_gpu<3,float, aggregate<float>>(phases.template get<0>(0).getDecomposition(),0);
    phases.add_no_device();
    phases.template get<0>(3) = vector_dist_gpu<3,float, aggregate<float>>(phases.template get<0>(0).getDecomposition(),0);
 
    for (size_t i = 0 ; i < 100 ; i++)
    {
        // Add a particle to all the phases
        phases.template get<0>(0).add();
        phases.template get<0>(1).add();
        phases.template get<0>(2).add();
        phases.template get<0>(3).add();
 
        phases.template get<0>(0).getLastPosWrite()[0] = (float)rand() / (float)RAND_MAX;
        phases.template get<0>(0).getLastPosWrite()[1] = (float)rand() / (float)RAND_MAX;
        phases.template get<0>(0).getLastPosWrite()[2] = (float)rand() / (float)RAND_MAX;
        phases.template get<0>(0).getLastPropWrite<0>() = (float)rand() / (float)RAND_MAX;
 
        phases.template get<0>(1).getLastPosWrite()[0] = (float)rand() / (float)RAND_MAX;
        phases.template get<0>(1).getLastPosWrite()[1] = (float)rand() / (float)RAND_MAX;
        phases.template get<0>(1).getLastPosWrite()[2] = (float)rand() / (float)RAND_MAX;
        phases.template get<0>(1).getLastPropWrite<0>() = (float)rand() / (float)RAND_MAX;
 
        phases.template get<0>(2).getLastPosWrite()[0] = (float)rand() / (float)RAND_MAX;
        phases.template get<0>(2).getLastPosWrite()[1] = (float)rand() / (float)RAND_MAX;
        phases.template get<0>(2).getLastPosWrite()[2] = (float)rand() / (float)RAND_MAX;
        phases.template get<0>(2).getLastPropWrite<0>() = (float)rand() / (float)RAND_MAX;
 
        phases.template get<0>(3).getLastPosWrite()[0] = (float)rand() / (float)RAND_MAX;
        phases.template get<0>(3).getLastPosWrite()[1] = (float)rand() / (float)RAND_MAX;
        phases.template get<0>(3).getLastPosWrite()[2] = (float)rand() / (float)RAND_MAX;
        phases.template get<0>(3).getLastPropWrite<0>() = (float)rand() / (float)RAND_MAX;
    }
 
    phases.template get<0>(0).hostToDevicePos();
    phases.template get<0>(1).hostToDevicePos();
    phases.template get<0>(2).hostToDevicePos();
    phases.template get<0>(3).hostToDevicePos();
 
    phases.template get<0>(0).template hostToDeviceProp<0>();
    phases.template get<0>(1).template hostToDeviceProp<0>();
    phases.template get<0>(2).template hostToDeviceProp<0>();
    phases.template get<0>(3).template hostToDeviceProp<0>();
 
    phases.template hostToDevice<0>();
 
    openfpm::vector_gpu<aggregate<float>> output;
    output.resize(100 * phases.size());
 
    CUDA_LAUNCH_DIM3(vdmkt,1,1,phases.toKernel(),output.toKernel());
 
    output.template deviceToHost<0>();
 
    bool match = true;
    int k = 0;
    for (int i = 0 ; i < phases.size() ; i++)
    {
        for (int j = 0 ; j < phases.template get<0>(i).size_local() ; j++)
        {
            match = output.template get<0>(k) == phases.template get<0>(i). template getProp<0>(j);
            k++;
        }
    }
 
    BOOST_REQUIRE_EQUAL(match,true);
}
 
BOOST_AUTO_TEST_CASE( vector_dist_multiphase_kernel_test_simplified )
{
    if (create_vcluster().getProcessingUnits() > 24)
    {return;}
 
    // The domain
    Box<3,float> box({-1000.0,-1000.0,-1000.0},{2000.0,2000.0,1000.0});
 
    // Boundary conditions
    size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
 
    float r_cut = 51.0;
 
    // ghost, big enough to contain the interaction radius
    Ghost<3,float> ghost(r_cut);
 
    openfpm::vector_custd<vector_dist_gpu<3,float, aggregate<float>> > phases;
 
    // first phase
    phases.add();
    phases.get(0) = vector_dist_gpu<3,float, aggregate<float>>(0,box,bc,ghost);
 
    // The other 3 phases
    phases.add();
    phases.get(1) = vector_dist_gpu<3,float, aggregate<float>>(phases.get(0).getDecomposition(),0);
    phases.add();
    phases.get(2) = vector_dist_gpu<3,float, aggregate<float>>(phases.get(0).getDecomposition(),0);
    phases.add();
    phases.get(3) = vector_dist_gpu<3,float, aggregate<float>>(phases.get(0).getDecomposition(),0);
 
    for (size_t i = 0 ; i < 100 ; i++)
    {
        // Add a particle to all the phases
        phases.get(0).add();
        phases.get(1).add();
        phases.get(2).add();
        phases.get(3).add();
 
        phases.get(0).getLastPosWrite()[0] = (float)rand() / (float)RAND_MAX;
        phases.get(0).getLastPosWrite()[1] = (float)rand() / (float)RAND_MAX;
        phases.get(0).getLastPosWrite()[2] = (float)rand() / (float)RAND_MAX;
        phases.get(0).getLastPropWrite<0>() = (float)rand() / (float)RAND_MAX;
 
        phases.get(1).getLastPosWrite()[0] = (float)rand() / (float)RAND_MAX;
        phases.get(1).getLastPosWrite()[1] = (float)rand() / (float)RAND_MAX;
        phases.get(1).getLastPosWrite()[2] = (float)rand() / (float)RAND_MAX;
        phases.get(1).getLastPropWrite<0>() = (float)rand() / (float)RAND_MAX;
 
        phases.get(2).getLastPosWrite()[0] = (float)rand() / (float)RAND_MAX;
        phases.get(2).getLastPosWrite()[1] = (float)rand() / (float)RAND_MAX;
        phases.get(2).getLastPosWrite()[2] = (float)rand() / (float)RAND_MAX;
        phases.get(2).getLastPropWrite<0>() = (float)rand() / (float)RAND_MAX;
 
        phases.get(3).getLastPosWrite()[0] = (float)rand() / (float)RAND_MAX;
        phases.get(3).getLastPosWrite()[1] = (float)rand() / (float)RAND_MAX;
        phases.get(3).getLastPosWrite()[2] = (float)rand() / (float)RAND_MAX;
        phases.get(3).getLastPropWrite<0>() = (float)rand() / (float)RAND_MAX;
    }
 
    phases.get(0).hostToDevicePos();
    phases.get(1).hostToDevicePos();
    phases.get(2).hostToDevicePos();
    phases.get(3).hostToDevicePos();
 
    phases.get(0).template hostToDeviceProp<0>();
    phases.get(1).template hostToDeviceProp<0>();
    phases.get(2).template hostToDeviceProp<0>();
    phases.get(3).template hostToDeviceProp<0>();
 
    phases.template hostToDevice<0>();
 
    openfpm::vector_gpu<aggregate<float>> output;
    output.resize(100 * phases.size());
 
    CUDA_LAUNCH_DIM3(vdmkt_simple,1,1,phases.toKernel(),output.toKernel());
 
    output.template deviceToHost<0>();
 
    bool match = true;
    int k = 0;
    for (int i = 0 ; i < phases.size() ; i++)
    {
        for (int j = 0 ; j < phases.get(i).size_local() ; j++)
        {
            match = output.template get<0>(k) == phases.get(i). template getProp<0>(j);
            k++;
        }
    }
 
    BOOST_REQUIRE_EQUAL(match,true);
}
 
BOOST_AUTO_TEST_CASE( vector_dist_multiphase_kernel_cl_test )
{
    if (create_vcluster().getProcessingUnits() > 24)
    {return;}
 
    // The domain
    Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
 
    // Boundary conditions
    size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
 
    // ghost, big enough to contain the interaction radius
    Ghost<3,float> ghost(0.1);
 
    openfpm::vector_custd<vector_dist_gpu<3,float, aggregate<float>> > phases;
    openfpm::vector_custd<decltype(phases.at(0).getCellListGPU(0.1))> cl_ph;
 
    // first phase
    phases.add();
    phases.get(0) = vector_dist_gpu<3,float, aggregate<float>>(0,box,bc,ghost);
 
    // The other 3 phases
    phases.add();
    phases.get(1) = vector_dist_gpu<3,float, aggregate<float>>(phases.get(0).getDecomposition(),0);
    phases.add();
    phases.get(2) = vector_dist_gpu<3,float, aggregate<float>>(phases.get(0).getDecomposition(),0);
    phases.add();
    phases.get(3) = vector_dist_gpu<3,float, aggregate<float>>(phases.get(0).getDecomposition(),0);
 
    for (size_t i = 0 ; i < 100 ; i++)
    {
        // Add a particle to all the phases
        phases.get(0).add();
        phases.get(1).add();
        phases.get(2).add();
        phases.get(3).add();
 
        phases.get(0).getLastPosWrite()[0] = (float)rand() / (float)RAND_MAX;
        phases.get(0).getLastPosWrite()[1] = (float)rand() / (float)RAND_MAX;
        phases.get(0).getLastPosWrite()[2] = (float)rand() / (float)RAND_MAX;
        phases.get(0).getLastPropWrite<0>() = (float)rand() / (float)RAND_MAX;
 
        phases.get(1).getLastPosWrite()[0] = (float)rand() / (float)RAND_MAX;
        phases.get(1).getLastPosWrite()[1] = (float)rand() / (float)RAND_MAX;
        phases.get(1).getLastPosWrite()[2] = (float)rand() / (float)RAND_MAX;
        phases.get(1).getLastPropWrite<0>() = (float)rand() / (float)RAND_MAX;
 
        phases.get(2).getLastPosWrite()[0] = (float)rand() / (float)RAND_MAX;
        phases.get(2).getLastPosWrite()[1] = (float)rand() / (float)RAND_MAX;
        phases.get(2).getLastPosWrite()[2] = (float)rand() / (float)RAND_MAX;
        phases.get(2).getLastPropWrite<0>() = (float)rand() / (float)RAND_MAX;
 
        phases.get(3).getLastPosWrite()[0] = (float)rand() / (float)RAND_MAX;
        phases.get(3).getLastPosWrite()[1] = (float)rand() / (float)RAND_MAX;
        phases.get(3).getLastPosWrite()[2] = (float)rand() / (float)RAND_MAX;
        phases.get(3).getLastPropWrite<0>() = (float)rand() / (float)RAND_MAX;
    }
 
    // redistribute all
 
    size_t tot = 0;
    size_t tot_g = 0;
    for (size_t i = 0 ; i < phases.size() ; i++)
    {
        phases.get(i).hostToDevicePos();
        phases.get(i).template hostToDeviceProp<0>();
        phases.get(i).map(RUN_ON_DEVICE);
        phases.get(i).template ghost_get<>(RUN_ON_DEVICE);
        phases.get(i).deviceToHostPos();
        phases.get(i).template deviceToHostProp<0>();
        tot += phases.get(i).size_local();
        tot_g += phases.get(i).size_local_with_ghost();
    }
 
    // Add cell list to the vector
 
    for (size_t i = 0 ; i < phases.size() ; i++)
    {
        cl_ph.add(phases.get(i).getCellListGPU(0.1));
    }
 
    //
 
    phases.template hostToDevice();
    cl_ph.template hostToDevice();
 
    openfpm::vector_gpu<aggregate<float>> output;
    openfpm::vector_gpu<aggregate<float>> output2;
    output.resize(tot);
    output2.resize(tot_g);
 
    CUDA_LAUNCH_DIM3(vdmkt_simple_cl,1,1,phases.toKernel(),output.toKernel(),cl_ph.toKernel(),output2.toKernel());
 
    output.template deviceToHost<0>();
 
    bool match = true;
    int k = 0;
    for (int i = 0 ; i < phases.size() ; i++)
    {
        for (int j = 0 ; j < phases.get(i).size_local() ; j++)
        {
            match = output.template get<0>(k) == phases.get(i). template getProp<0>(j);
            k++;
        }
    }
 
    BOOST_REQUIRE_EQUAL(match,true);
}
 
BOOST_AUTO_TEST_SUITE_END()