Commit d8f45e69 authored by incardon's avatar incardon

Fixing Cell-list for GPU + tests + other small fixes

parent 7ae9fb58
......@@ -49,7 +49,7 @@ struct copy_cpu_encap_encap_prp
__device__ __host__ inline copy_cpu_encap_encap_prp(e_src & src, e_dst & dst)
:src(src),dst(dst)
{
#ifdef SE_CLASS1
#if defined(SE_CLASS1) && !defined(__NVCC__)
// e_src and e_dst must have the same number of properties
if (e_src::T_type::max_prop != e_dst::T_type::max_prop)
......@@ -111,7 +111,133 @@ struct copy_cpu_encap_encap
__device__ __host__ inline copy_cpu_encap_encap(const e_src & src, e_dst & dst)
:src(src),dst(dst)
{
#if defined(SE_CLASS1) && !defined(__NVCC__)
// e_src and e_dst must have the same number of properties
if (e_src::T_type::max_prop != e_dst::T_type::max_prop)
std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " the number of properties between src and dst must match";
#endif
};
#ifdef SE_CLASS1
/*! \brief Constructor
*
* Calling this constructor produce an error. This class store the reference of the object,
* this mean that the object passed must not be a temporal object
*
*/
inline copy_cpu_encap_encap(const e_src && src, const e_dst && dst)
:src(src),dst(dst)
{std::cerr << "Error: " <<__FILE__ << ":" << __LINE__ << " Passing a temporal object";};
#endif
//! It call the copy function for each property
template<typename T>
__device__ __host__ inline void operator()(T& t) const
{
// Remove the reference from the type to copy
typedef typename boost::remove_reference<decltype(dst.template get<T::value>())>::type copy_rtype;
meta_copy<copy_rtype>::meta_copy_(src.template get<T::value>(),dst.template get<T::value>());
}
};
/*! \brief this class is a functor for "for_each" algorithm
*
* This class is a functor for "for_each" algorithm. For each
* element of the boost::vector the operator() is called.
* Is mainly used to copy one encap into another encap object
*
* \tparam encap source
* \tparam encap dst
*
*/
template<typename e_src, typename e_dst>
struct copy_cpu_encap_encap_general
{
//! encapsulated source object
const e_src & src;
//! encapsulated destination object
e_dst & dst;
/*! \brief constructor
*
* \param src source encapsulated object
* \param dst source encapsulated object
*
*/
__device__ __host__ inline copy_cpu_encap_encap_general(const e_src & src, e_dst & dst)
:src(src),dst(dst)
{
#if defined(SE_CLASS1) && !defined(__NVCC__)
// e_src and e_dst must have the same number of properties
if (e_src::T_type::max_prop != e_dst::T_type::max_prop)
std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " the number of properties between src and dst must match";
#endif
};
#ifdef SE_CLASS1
/*! \brief Constructor
*
* Calling this constructor produce an error. This class store the reference of the object,
* this mean that the object passed must not be a temporal object
*
*/
inline copy_cpu_encap_encap_general(const e_src && src, const e_dst && dst)
:src(src),dst(dst)
{std::cerr << "Error: " <<__FILE__ << ":" << __LINE__ << " Passing a temporal object";};
#endif
//! It call the copy function for each property
template<typename T>
__device__ __host__ inline void operator()(T& t) const
{
// Remove the reference from the type to copy
typedef typename boost::remove_reference<decltype(dst.template get<T::value>())>::type copy_dtype;
// Remove the reference from the type to copy
typedef typename boost::remove_reference<decltype(src.template get<T::value>())>::type copy_stype;
meta_copy_d<copy_stype,copy_dtype>::meta_copy_d_(src.template get<T::value>(),dst.template get<T::value>());
}
};
/*! \brief this class is a functor for "for_each" algorithm
*
* This class is a functor for "for_each" algorithm. For each
* element of the boost::vector the operator() is called.
* Is mainly used to copy one encap into another encap object
*
* \tparam encap source
* \tparam encap dst
*
*/
template<typename e_src>
struct copy_cpu_encap_single
{
//! encapsulated source object
const e_src & src;
//! encapsulated destination object
e_src & dst;
/*! \brief constructor
*
* \param src source encapsulated object
* \param dst source encapsulated object
*
*/
__device__ __host__ inline copy_cpu_encap_single(const e_src & src, e_src & dst)
:src(src),dst(dst)
{
#if defined(SE_CLASS1) && !defined(__NVCC__)
// e_src and e_dst must have the same number of properties
if (e_src::T_type::max_prop != e_dst::T_type::max_prop)
......@@ -143,6 +269,7 @@ struct copy_cpu_encap_encap
}
};
///////////////////////////////////
//! It copy two encap object
......@@ -425,7 +552,7 @@ public:
*/
__device__ __host__ inline encapc<dim,T,Mem> & operator=(const encapc<dim,T,Mem> & ec)
{
copy_cpu_encap_encap<encapc<dim,T,Mem>,encapc<dim,T,Mem>> cp(ec,*this);
copy_cpu_encap_single<encapc<dim,T,Mem>> cp(ec,*this);
boost::mpl::for_each_ref< boost::mpl::range_c<int,0,T::max_prop> >(cp);
......@@ -441,7 +568,7 @@ public:
*/
__device__ __host__ inline encapc<dim,T,Mem> & operator=(const encapc<dim,T,Mem2> & ec)
{
copy_cpu_encap_encap<encapc<dim,T,Mem2>,encapc<dim,T,Mem>> cp(ec,*this);
copy_cpu_encap_encap_general<encapc<dim,T,Mem2>,encapc<dim,T,Mem>> cp(ec,*this);
boost::mpl::for_each_ref< boost::mpl::range_c<int,0,T::max_prop> >(cp);
......@@ -599,7 +726,7 @@ public:
*/
__device__ __host__ inline encapc<dim,T,Mem> & operator=(const encapc<dim,T,Mem> & ec)
{
copy_cpu_encap_encap<encapc<dim,T,Mem>,encapc<dim,T,Mem>> cp(ec,*this);
copy_cpu_encap_single<encapc<dim,T,Mem>> cp(ec,*this);
boost::mpl::for_each_ref< boost::mpl::range_c<int,0,T::max_prop> >(cp);
......@@ -615,7 +742,7 @@ public:
*/
__device__ __host__ inline encapc<dim,T,Mem> & operator=(const encapc<dim,T,Mem2> & ec)
{
copy_cpu_encap_encap<encapc<dim,T,Mem2>,encapc<dim,T,Mem>> cp(ec,*this);
copy_cpu_encap_encap_general<encapc<dim,T,Mem2>,encapc<dim,T,Mem>> cp(ec,*this);
boost::mpl::for_each_ref< boost::mpl::range_c<int,0,T::max_prop> >(cp);
......
......@@ -100,15 +100,15 @@ struct grid_gpu_ker
//! layout data
layout data_;
grid_gpu_ker()
__device__ __host__ grid_gpu_ker()
{}
grid_gpu_ker(const grid_sm<dim,void> & g1)
__device__ __host__ grid_gpu_ker(const grid_sm<dim,void> & g1)
:g1(g1)
{
}
grid_gpu_ker(const grid_gpu_ker & cpy)
__device__ __host__ grid_gpu_ker(const grid_gpu_ker & cpy)
:g1(cpy.g1)
{
grid_gpu_ker_constructor_impl<is_layout_inte<layout_base<T_>>::value,T_>::construct(cpy,*this);
......
......@@ -241,11 +241,13 @@ private:
*/
inline void check_init() const
{
#ifndef __NVCC__
if (is_mem_init == false)
{
std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " you must call SetMemory before access the grid\n";
ACTION_ON_ERROR(GRID_ERROR_OBJECT);
}
#endif
}
/*! \brief Check that the key is inside the grid
......@@ -255,6 +257,7 @@ private:
*/
inline void check_bound(const grid_key_dx<dim> & v1) const
{
#ifndef __NVCC__
for (long int i = 0 ; i < dim ; i++)
{
if (v1.get(i) >= (long int)getGrid().size(i))
......@@ -268,6 +271,7 @@ private:
ACTION_ON_ERROR(GRID_ERROR_OBJECT);
}
}
#endif
}
/*! \brief Check that the key is inside the grid
......@@ -277,11 +281,13 @@ private:
*/
inline void check_bound(size_t v1) const
{
#ifndef __NVCC__
if (v1 >= getGrid().size())
{
std::cerr << "Error " __FILE__ << ":" << __LINE__ <<" grid overflow " << v1<< " >= " << getGrid().size() << "\n";
ACTION_ON_ERROR(GRID_ERROR_OBJECT);
}
#endif
}
/*! \brief Check that the key is inside the grid
......@@ -294,6 +300,35 @@ private:
*/
template<typename Mem> inline void check_bound(const grid_base_impl<dim,T,Mem,layout,layout_base> & g,const grid_key_dx<dim> & key2) const
{
#ifndef __NVCC__
for (size_t i = 0 ; i < dim ; i++)
{
if (key2.get(i) >= (long int)g.getGrid().size(i))
{
std::cerr << "Error " __FILE__ << ":" << __LINE__ <<" grid overflow " << "x=[" << i << "]=" << key2.get(i) << " >= " << g.getGrid().size(i) << "\n";
ACTION_ON_ERROR(GRID_ERROR_OBJECT);
}
else if (key2.get(i) < 0)
{
std::cerr << "Error " __FILE__ << ":" << __LINE__ <<" grid overflow " << "x=[" << i << "]=" << key2.get(i) << " is negative " << "\n";
ACTION_ON_ERROR(GRID_ERROR_OBJECT);
}
}
#endif
}
/*! \brief Check that the key is inside the grid
*
* check if key2 is inside the g grid boundary
*
* \param g grid
* \param key2
*
*/
template<typename Mem, typename layout2, template <typename>
class layout_base2> inline void check_bound(const grid_base_impl<dim,T,Mem,layout2,layout_base2> & g,const grid_key_dx<dim> & key2) const
{
#ifndef __NVCC__
for (size_t i = 0 ; i < dim ; i++)
{
if (key2.get(i) >= (long int)g.getGrid().size(i))
......@@ -307,6 +342,7 @@ private:
ACTION_ON_ERROR(GRID_ERROR_OBJECT);
}
}
#endif
}
/*! \brief Check that the key is inside the grid
......@@ -319,11 +355,13 @@ private:
*/
template<typename Mem> inline void check_bound(const grid_base_impl<dim,T,Mem,layout,layout_base> & g,const size_t & key2) const
{
#ifndef __NVCC__
if (key2 >= g.getGrid().size())
{
std::cerr << "Error " __FILE__ << ":" << __LINE__ <<" grid overflow " << key2 << " >= " << getGrid().size() << "\n";
ACTION_ON_ERROR(GRID_ERROR_OBJECT);
}
#endif
}
#endif
......
......@@ -102,7 +102,7 @@ public:
*/
template<typename ...T> __device__ __host__ inline grid_key_dx(const size_t v,const T...t)
{
#ifdef SE_CLASS1
#if defined(SE_CLASS1) && !defined(__NVCC__)
if (sizeof...(t) != dim -1)
{std::cerr << "Error grid_key: " << __FILE__ << " " << __LINE__ << " creating a key of dimension " << dim << " require " << dim << " numbers " << sizeof...(t) + 1 << " provided" << "\n";}
#endif
......@@ -415,7 +415,7 @@ public:
*/
__device__ __host__ void set_d(size_t i, mem_id id)
{
#ifdef SE_CLASS1
#if defined(SE_CLASS1) && !defined(__NVCC__)
if (i >= dim)
std::cerr << "grid_key_dx error: " << __FILE__ << " " << __LINE__ << " try to access dimension " << i << " on a grid_key_dx of size " << dim << "\n";
......
......@@ -212,6 +212,90 @@ template<unsigned int dim> void NNcalc_full(openfpm::vector<grid_key_dx<dim>> &
};
/*! Calculate the neighborhood cells based on the radius
*
* \note To the calculated neighborhood cell you have to add the id of the central cell
*
\verbatim
+-----------------------+
|p |p |p |p |p |p |p |p |
+-----------------------+
|p | | | | | | |p |
+-----------------------+
|p | | |7 |8 |9 | |p |
+-----------------------+
|p | | |-1|0 |1 | |p |
+-----------------------+
|p |9 | |-9|-8|-7| |p |
+-----------------------+
|p |p |p |p |p |p |p |p |
+-----------------------+
\endverbatim
*
* The number indicate the cell id calculated
*
* -9,-8,-7,-1,0,1,7,8,9
*
* The cell 0 has id = 22 in the big cell matrix, so to calculate the
* neighborhood cells you have to sum the id of the center cell
*
* 13,14,15,21,22,23,29,30,31
*
* \param r_cut Cutoff-radius
* \param NNcell vector containing the neighborhood cells ids
*
*/
template<unsigned int dim, typename T>
void NNcalc_rad(T r_cut, openfpm::vector<long int> & NNcell, const Box<dim,T> & cell_box, const grid_sm<dim,void> & gs)
{
size_t n_cell[dim];
size_t n_cell_mid[dim];
Point<dim,T> spacing = cell_box.getP2();
for (size_t i = 0 ; i < dim ; i++)
{
n_cell[i] = 2*(std::ceil(r_cut / spacing.get(i)))+1;
n_cell_mid[i] = n_cell[i] / 2;
}
grid_sm<dim,void> gsc(n_cell);
grid_key_dx_iterator<dim> gkdi(gsc);
Box<dim,T> cell_zero;
for (unsigned int i = 0 ; i < dim ; i++)
{
cell_zero.setLow(i,n_cell_mid[i]*spacing.get(i));
cell_zero.setHigh(i,(n_cell_mid[i]+1)*spacing.get(i));
}
while (gkdi.isNext())
{
auto key = gkdi.get();
Box<dim,T> cell;
for (unsigned int i = 0 ; i < dim ; i++)
{
cell.setLow(i,key.get(i)*spacing.get(i));
cell.setHigh(i,(key.get(i)+1)*spacing.get(i));
}
// here we check if the cell is in the radius.
T min_distance = cell.min_distance(cell_zero);
if (min_distance > r_cut)
{++gkdi;continue;}
for (size_t i = 0 ; i < dim ; i++)
{key.set_d(i,key.get(i) - n_cell_mid[i]);}
NNcell.add(gs.LinId(key));
++gkdi;
}
}
/* NOTE all the implementations
*
* has complexity O(1) in getting the cell id and the elements in a cell
......@@ -300,68 +384,10 @@ private:
//! has been constructed from an old decomposition
size_t n_dec;
/*! Calculate the neighborhood cells based on the radius
*
* \note To the calculated neighborhood cell you have to add the id of the central cell
*
\verbatim
+-----------------------+
|p |p |p |p |p |p |p |p |
+-----------------------+
|p | | | | | | |p |
+-----------------------+
|p | | |7 |8 |9 | |p |
+-----------------------+
|p | | |-1|0 |1 | |p |
+-----------------------+
|p |9 | |-9|-8|-7| |p |
+-----------------------+
|p |p |p |p |p |p |p |p |
+-----------------------+
\endverbatim
*
* The number indicate the cell id calculated
*
* -9,-8,-7,-1,0,1,7,8,9
*
* The cell 0 has id = 22 in the big cell matrix, so to calculate the
* neighborhood cells you have to sum the id of the center cell
*
* 13,14,15,21,22,23,29,30,31
*
* \param r_cut Cutoff-radius
* \param NNcell vector containing the neighborhood cells ids
*
*/
void NNcalc(T r_cut, openfpm::vector<long int> & NNcell)
{
size_t n_cell[dim];
size_t n_cell_mid[dim];
Point<dim,T> spacing = this->getCellBox().getP2();
const grid_sm<dim,void> & gs = this->getGrid();
//! Cells for the neighborhood radius
openfpm::vector<long int> nnc_rad;
for (size_t i = 0 ; i < dim ; i++)
{
n_cell[i] = 2*(std::ceil(r_cut / spacing.get(i)))+1;
n_cell_mid[i] = n_cell[i] / 2;
}
grid_sm<dim,void> gsc(n_cell);
grid_key_dx_iterator<dim> gkdi(gsc);
while (gkdi.isNext())
{
auto key = gkdi.get();
for (size_t i = 0 ; i < dim ; i++)
key.set_d(i,key.get(i) - n_cell_mid[i]);
NNcell.add(gs.LinId(key));
++gkdi;
}
}
//! Initialize the structures of the data structure
void InitializeStructures(const size_t (& div)[dim], size_t tot_n_cell, size_t slot=STARTING_NSLOT)
......@@ -652,6 +678,16 @@ public:
return *this;
}
/*! \brief Set the radius for the getNNIteratorRadius
*
* \param radius
*
*/
void setRadius(T radius)
{
NNcalc_rad(radius,nnc_rad,this->getCellBox(),this->getGrid());
}
/*! \brief Get an iterator over particles following the cell structure
*
* \param dom_cells cells in the domain
......@@ -793,6 +829,15 @@ public:
Mem_type::remove(cell,ele);
}
/*! \brief Get the number of cells this cell-list contain
*
* \return number of cells
*/
inline size_t getNCells() const
{
return Mem_type::size();
}
/*! \brief Return the number of elements in the cell
*
* \param cell_id id of the cell
......@@ -862,6 +907,32 @@ public:
return CellIterator<CellList<dim,T,Mem_type,transform>>(cell,*this);
}
/*! \brief Get the Neighborhood iterator
*
* It iterate across all the element of the selected cell and the near cells
*
* \verbatim
* * *
* x *
* * *
\endverbatim
*
* * x is the selected cell
* * * are the near cell
*
* \param cell cell id
*
* \return An iterator across the neighhood particles
*
*/
template<unsigned int impl=NO_CHECK> inline CellNNIteratorRadius<dim,CellList<dim,T,Mem_type,transform>,impl> getNNIteratorRadius(size_t cell)
{
CellNNIteratorRadius<dim,CellList<dim,T,Mem_type,transform>,impl> cln(cell,nnc_rad,*this);
return cln;
}
/*! \brief Get the Neighborhood iterator
*
* It iterate across all the element of the selected cell and the near cells
......@@ -904,7 +975,7 @@ public:
openfpm::vector<long int> & NNc = rcache[r_cut];
if (NNc.size() == 0)
NNcalc(r_cut,NNc);
{NNcalc_rad(r_cut,NNc,this->getCellBox(),this->getGrid());}
CellNNIteratorRadius<dim,CellList<dim,T,Mem_type,transform>,impl> cln(cell,NNc,*this);
......
This diff is collapsed.
......@@ -351,8 +351,160 @@ template<typename CellList> void Test_CellDecomposer_consistent()
BOOST_REQUIRE_EQUAL(bx_sub.getHigh(1),1.0f/5.0f);
}
/*! \brief Test cell structure
*
* \tparam CellS
*
*/
template<unsigned int dim, typename T, typename CellS> void Test_NN_iterator_radius(SpaceBox<dim,T> & box)
{
//Space where is living the Cell list
//SpaceBox<dim,T> box({0.0f,0.0f,0.0f},{1.0f,1.0f,1.0f});
// Subdivisions
size_t div[dim];
size_t div2[dim];
for (size_t i = 0 ; i < dim ; i++)
{
div[i] = 17;
div2[i] = 34;
}
// grid info
grid_sm<dim,void> g_info(div);
grid_sm<dim,void> g_info2(div2);
//! [Usage of cell list multi]
// CellS = CellListM<dim,T,8>
CellS cl1(box,div);
CellS cl2(box,div2,2);
T radius = (box.getHigh(0) - box.getLow(0))/div[0];
cl2.setRadius( radius );
// create a vector of random point
openfpm::vector<Point<dim,float>> vrp;
for (size_t j = 0 ; j < 10000 ; j++)
{
vrp.add();
for (size_t i = 0 ; i < dim ; i++)
{
vrp.template get<0>(j)[i] = ((float)rand() / RAND_MAX)*(box.getHigh(i) - box.getLow(i)) + box.getLow(i);
}
}
auto g_it = vrp.getIterator();
while (g_it.isNext())
{
Point<dim,T> xp = vrp.get(g_it.get());
size_t debug = cl1.getCell(xp);
cl1.add(xp,g_it.get());
cl2.add(xp,g_it.get());
++g_it;
}
// Get the neighborhood of each particle and compare the numbers of cell
bool match = true;
auto g_it2 = vrp.getIterator();
size_t number_of_nn = 0;
size_t number_of_nn2 = 0;