diff --git a/CHANGELOG.md b/CHANGELOG.md
index 28586c897cea1be9c9d26e64910638286ffe11ef..fa05c4e5467b391523890a1e7286a657479402b9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,7 +1,19 @@
# Change Log
All notable changes to this project will be documented in this file.
-## [0.5.0 - Gingold] - Mid August 2016
+## [0.5.1] - End of August
+
+### Added
+- Symmetric cell list
+- Verlet list
+- Full-Support for complex property on vector-dist (Serialization) + example
+
+### Fixed
+- Energy calculation in md example (double counting potential energy)
+
+### Changed
+
+## [0.5.0] - 15 August 2016
### Added
- map communicate particles across processors mooving the information of all the particle map_list give the possibility to give a list of property to move from one to another processor
diff --git a/example/Vector/3_molecular_dynamic/main.cpp b/example/Vector/3_molecular_dynamic/main.cpp
index ac62885fbfb32f193be2b4206fac2222888b6bea..81188e97d753597ca1bc9cfeec51cb8fbd127e85 100644
--- a/example/Vector/3_molecular_dynamic/main.cpp
+++ b/example/Vector/3_molecular_dynamic/main.cpp
@@ -218,11 +218,6 @@ double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd,
// Get the position of the particle p
Point<3,double> xp = vd.getPos(p);
- // Reset the force
- vd.template getProp<force>(p)[0] = 0.0;
- vd.template getProp<force>(p)[1] = 0.0;
- vd.template getProp<force>(p)[2] = 0.0;
-
// Get an iterator over the neighborhood of the particle p
auto Np = NN.template getNNIterator<NO_CHECK>(NN.getCell(vd.getPos(p)));
@@ -242,7 +237,7 @@ double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd,
double rn = norm2(xp - xq);
// potential energy (using pow is slower)
- E += 4.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) );
+ E += 2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) );
// Next neighborhood
++Np;
@@ -286,8 +281,8 @@ int main(int argc, char* argv[])
//! \cond [constants] \endcond
double dt = 0.0005;
- float r_cut = 0.3;
double sigma = 0.1;
+ double r_cut = 3.0*sigma;
double sigma12 = pow(sigma,12);
double sigma6 = pow(sigma,6);
diff --git a/example/Vector/3_molecular_dynamic/main_expr.cpp b/example/Vector/3_molecular_dynamic/main_expr.cpp
index 86c7e31bf41bbd94ef81c469fa90617ed261b0b1..57f474efe8f2c12ffd35f2f253e8d4afab28f348 100644
--- a/example/Vector/3_molecular_dynamic/main_expr.cpp
+++ b/example/Vector/3_molecular_dynamic/main_expr.cpp
@@ -15,7 +15,7 @@ struct ln_potential
{
double rn = norm2(xp - xq);
- Point<2,double> E({4.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ),
+ Point<2,double> E({2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ),
0.0});
return E;
diff --git a/openfpm_vcluster b/openfpm_vcluster
index a874640c9b69979cbf3117a4ecf9e471ae435b62..0039a8c1bc37029fdf99aaf633b1dfa64071fffc 160000
--- a/openfpm_vcluster
+++ b/openfpm_vcluster
@@ -1 +1 @@
-Subproject commit a874640c9b69979cbf3117a4ecf9e471ae435b62
+Subproject commit 0039a8c1bc37029fdf99aaf633b1dfa64071fffc
diff --git a/src/Grid/grid_dist_id.hpp b/src/Grid/grid_dist_id.hpp
index 5450c9e61ce8b21e20192e2df44bc798278da6a1..c16ae0f8fdf5843aabdeee22949c42995a31b2a1 100644
--- a/src/Grid/grid_dist_id.hpp
+++ b/src/Grid/grid_dist_id.hpp
@@ -519,15 +519,9 @@ class grid_dist_id
// check that the grid has valid size
check_size(g_sz);
- // For a 5x5 grid you have 4x4 Cell (With the exception of periodic)
+ // get the size of the cell decomposer
size_t c_g[dim];
- for (size_t i = 0 ; i < dim ; i++)
- {
- if (bc[i] == NON_PERIODIC)
- c_g[i] = (g_sz[i]-1 > 0)?(g_sz[i]-1):1;
- else
- c_g[i] = g_sz[i];
- }
+ getCellDecomposerPar<dim>(c_g,g_sz,bc);
// Initialize the cell decomposer
cd_sm.setDimensions(domain,c_g,0);
diff --git a/src/Grid/grid_dist_id_iterator_dec.hpp b/src/Grid/grid_dist_id_iterator_dec.hpp
index f63addc735d17b26b73f6d69b5a493e1a6481f54..c9425784d42be6817b2a40dcaa97054418ea8d23 100644
--- a/src/Grid/grid_dist_id_iterator_dec.hpp
+++ b/src/Grid/grid_dist_id_iterator_dec.hpp
@@ -145,6 +145,7 @@ class grid_dist_id_iterator_dec
*
* \param dec Decomposition
* \param sz size of the grid
+ * \param bc boundary conditions
*
*/
grid_dist_id_iterator_dec(Decomposition & dec, const size_t (& sz)[Decomposition::dims])
@@ -152,7 +153,8 @@ class grid_dist_id_iterator_dec
{
// Initialize start and stop
start.zero();
- for (size_t i = 0 ; i < Decomposition::dims ; i++) stop.set_d(i,sz[i]-1);
+ for (size_t i = 0 ; i < Decomposition::dims ; i++)
+ stop.set_d(i,sz[i]-1);
// From the decomposition construct gdb_ext
create_gdb_ext<Decomposition::dims,Decomposition>(gdb_ext,dec,sz,dec.getDomain(),spacing);
@@ -211,27 +213,6 @@ class grid_dist_id_iterator_dec
return *this;
}
- /*! \brief Convert a g_dist_key_dx into a global key
- *
- * \see grid_dist_key_dx
- * \see grid_dist_iterator
- *
- * \return the global position in the grid
- *
- */
-/* inline grid_key_dx<Decomposition::dims> getGKey(const grid_dist_key_dx<Decomposition::dims> & k)
- {
- // Get the sub-domain id
- size_t sub_id = k.getSub();
-
- grid_key_dx<Decomposition::dims> k_glob = k.getKey();
-
- // shift
- k_glob = k_glob + gdb_ext.get(sub_id).origin;
-
- return k_glob;
- }*/
-
/*! \brief Check if there is the next element
*
* \return true if there is the next, false otherwise
diff --git a/src/Grid/grid_dist_util.hpp b/src/Grid/grid_dist_util.hpp
index 2981c19900f97e5bf192c353fcf0d003faf5bcee..bba37eb0335b80bd7e99aeb03f9be68a185861d4 100644
--- a/src/Grid/grid_dist_util.hpp
+++ b/src/Grid/grid_dist_util.hpp
@@ -10,6 +10,25 @@
#include "NN/CellList/CellDecomposer.hpp"
+/*! \brief get cellDecomposer parameters
+ *
+ * \tparam dim dimensionality
+ *
+ * \param c_g get the parameters of the cell decomposer
+ * \param g_sz global grid parameters
+ *
+ */
+template<unsigned int dim> void getCellDecomposerPar(size_t (& c_g)[dim], const size_t (& g_sz)[dim], const size_t (& bc)[dim])
+{
+ for (size_t i = 0 ; i < dim ; i++)
+ {
+ if (bc[i] == NON_PERIODIC)
+ c_g[i] = (g_sz[i]-1 > 0)?(g_sz[i]-1):1;
+ else
+ c_g[i] = g_sz[i];
+ }
+}
+
/*! \brief
*
*
@@ -87,20 +106,21 @@ template<int dim, typename Decomposition> inline void create_gdb_ext(openfpm::ve
* \param sz Global grid grid size
* \param domain Domain where the grid is defined
* \param spacing Define the spacing of the grid
+ * \param bc boundary conditions
*
*/
template<int dim, typename Decomposition> inline void create_gdb_ext(openfpm::vector<GBoxes<dim>> & gdb_ext, Decomposition & dec, const size_t (& sz)[dim], const Box<Decomposition::dims,typename Decomposition::stype> & domain, typename Decomposition::stype (& spacing)[dim])
{
// Create the cell decomposer
-
CellDecomposer_sm<Decomposition::dims,typename Decomposition::stype, shift<Decomposition::dims,typename Decomposition::stype>> cd_sm;
- size_t sz_cell[Decomposition::dims];
- for (size_t i = 0 ; i < dim ; i++)
- sz_cell[i] = sz[i] - 1;
+ size_t cdp[dim];
+
+ // Get the parameters to create a Cell-decomposer
+ getCellDecomposerPar<Decomposition::dims>(cdp,sz,dec.periodicity());
// Careful cd_sm require the number of cell
- cd_sm.setDimensions(domain,sz_cell,0);
+ cd_sm.setDimensions(domain,cdp,0);
create_gdb_ext<dim,Decomposition>(gdb_ext,dec,cd_sm);
diff --git a/src/Vector/vector_dist.hpp b/src/Vector/vector_dist.hpp
index 309563ee13a466d7e05b4c1b1489bf0489aeb0b9..f5348417c3882d42b06ff68b1443fae8d1017143 100644
--- a/src/Vector/vector_dist.hpp
+++ b/src/Vector/vector_dist.hpp
@@ -29,6 +29,7 @@
#include "Vector/vector_dist_ofb.hpp"
#include "Decomposition/CartDecomposition.hpp"
#include "data_type/aggregate.hpp"
+#include "NN/VerletList/VerletList.hpp"
#define V_SUB_UNIT_FACTOR 64
@@ -776,8 +777,7 @@ private:
*
* \return the processor bounding box
*/
- inline Box<dim, St> cl_param_calculate(size_t (&div)[dim], St r_cut, const Ghost<dim, St> & enlarge)
-
+/* inline Box<dim, St> cl_param_calculate(size_t (&div)[dim], St r_cut, const Ghost<dim, St> & enlarge)
{
// calculate the parameters of the cell list
@@ -795,7 +795,7 @@ private:
pbox.setHigh(i,pbox.getLow(i) + div[i]*r_cut);
}
return pbox;
- }
+ }*/
/*! \brief Initialize the structures
*
@@ -1456,8 +1456,11 @@ public:
// Division array
size_t div[dim];
+ // get the processor bounding box
+ Box<dim, St> pbox = dec.getProcessorBounds();
+
// Processor bounding box
- auto pbox = cl_param_calculate(div, r_cut, enlarge);
+ cl_param_calculate(pbox, div, r_cut, enlarge);
cell_list.Initialize(pbox, div);
@@ -1488,8 +1491,11 @@ public:
// Division array
size_t div[dim];
+ // get the processor bounding box
+ Box<dim, St> pbox = dec.getProcessorBounds();
+
// Processor bounding box
- auto pbox = cl_param_calculate(div, r_cut, enlarge);
+ cl_param_calculate(pbox,div, r_cut, enlarge);
cell_list.Initialize(pbox, div, g_m);
@@ -1498,13 +1504,39 @@ public:
return cell_list;
}
+
+ /*! \brief for each particle get the verlet list
+ *
+ * \param r_cut cut-off radius
+ *
+ */
+ VerletList<dim,St,FAST,shift<dim,St> > getVerlet(St r_cut)
+ {
+ VerletList<dim,St,FAST,shift<dim,St>> ver;
+
+ // Division array
+ size_t div[dim];
+
+ // get the processor bounding box
+ Box<dim, St> bt = dec.getProcessorBounds();
+
+ // Calculate the divisions for the Cell-lists
+ cl_param_calculate(bt,div,r_cut,Ghost<dim,St>(0.0));
+
+ ver.Initialize(bt,r_cut,v_pos,g_m);
+
+ return ver;
+ }
+
/*! \brief for each particle get the verlet list
*
* \param verlet output verlet list for each particle
* \param r_cut cut-off radius
*
+ * \deprecated
+ *
*/
- void getVerlet(openfpm::vector<openfpm::vector<size_t>> & verlet, St r_cut)
+ void getVerletDeprecated(openfpm::vector<openfpm::vector<size_t>> & verlet, St r_cut)
{
// resize verlet to store the number of particles
verlet.resize(size_local());
@@ -1722,25 +1754,15 @@ public:
*/
inline grid_dist_id_iterator_dec<Decomposition> getGridIterator(const size_t (&sz)[dim])
{
- size_t sz_g[dim];
grid_key_dx<dim> start;
grid_key_dx<dim> stop;
for (size_t i = 0; i < dim; i++)
{
start.set_d(i, 0);
- if (dec.periodicity(i) == PERIODIC)
- {
- sz_g[i] = sz[i];
- stop.set_d(i, sz_g[i] - 2);
- }
- else
- {
- sz_g[i] = sz[i];
- stop.set_d(i, sz_g[i] - 1);
- }
+ stop.set_d(i, sz[i] - 1);
}
- grid_dist_id_iterator_dec<Decomposition> it_dec(dec, sz_g, start, stop);
+ grid_dist_id_iterator_dec<Decomposition> it_dec(dec, sz, start, stop);
return it_dec;
}
diff --git a/src/Vector/vector_dist_unit_test.hpp b/src/Vector/vector_dist_unit_test.hpp
index d5c12a1df3de34e19e6af10ded343934b36fc3dc..6c300d2643b599a2b4d47dedccf9136e41b549a6 100644
--- a/src/Vector/vector_dist_unit_test.hpp
+++ b/src/Vector/vector_dist_unit_test.hpp
@@ -1086,6 +1086,71 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_interacting_particles )
}
}
+BOOST_AUTO_TEST_CASE( vector_dist_grid_iterator )
+{
+ long int k = 64*64*64*create_vcluster().getProcessingUnits();
+ k = std::pow(k, 1/3.);
+
+ long int big_step = k / 30;
+ big_step = (big_step == 0)?1:big_step;
+ long int small_step = 21;
+
+ print_test( "Testing vector grid iterator list k<=",k);
+
+ // 3D test
+ for ( ; k > 8*big_step ; k-= (k > 2*big_step)?big_step:small_step )
+ {
+ Vcluster & v_cl = create_vcluster();
+
+ const size_t Ng = k;
+
+ // we create a 128x128x128 Grid iterator
+ size_t sz[3] = {Ng,Ng,Ng};
+
+ Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
+
+ // Boundary conditions
+ size_t bc[3]={NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
+
+ // ghost
+ Ghost<3,float> ghost(1.0/(Ng-2));
+
+ // Distributed vector
+ vector_dist<3,float, Point_test<float>, CartDecomposition<3,float> > vd(0,box,bc,ghost);
+
+ // Put particles on a grid creating a Grid iterator
+ auto it = vd.getGridIterator(sz);
+
+ while (it.isNext())
+ {
+ vd.add();
+
+ auto key = it.get();
+
+ vd.getLastPos()[0] = key.get(0) * it.getSpacing(0);
+ vd.getLastPos()[1] = key.get(1) * it.getSpacing(1);
+ vd.getLastPos()[2] = key.get(2) * it.getSpacing(2);
+
+ ++it;
+ }
+
+ BOOST_REQUIRE_EQUAL(it.getSpacing(0),1.0f/(Ng-1));
+ BOOST_REQUIRE_EQUAL(it.getSpacing(1),1.0f/(Ng-1));
+ BOOST_REQUIRE_EQUAL(it.getSpacing(2),1.0f/(Ng-1));
+
+ // distribute particles and sync ghost
+ vd.map();
+
+
+ // Check that the sum of all the particles is the grid size
+ size_t total = vd.size_local();
+ v_cl.sum(total);
+ v_cl.execute();
+
+ BOOST_REQUIRE_EQUAL(total,(Ng) * (Ng) * (Ng));
+ }
+}
+
BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
{
long int k = 64*64*64*create_vcluster().getProcessingUnits();
@@ -1134,6 +1199,10 @@ BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
++it;
}
+ BOOST_REQUIRE_EQUAL(it.getSpacing(0),1.0f/Ng);
+ BOOST_REQUIRE_EQUAL(it.getSpacing(1),1.0f/Ng);
+ BOOST_REQUIRE_EQUAL(it.getSpacing(2),1.0f/Ng);
+
// distribute particles and sync ghost
vd.map();
@@ -1142,7 +1211,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
v_cl.sum(total);
v_cl.execute();
- BOOST_REQUIRE_EQUAL(total,(Ng-1) * (Ng-1) * (Ng-1));
+ BOOST_REQUIRE_EQUAL(total,(Ng) * (Ng) * (Ng));
vd.ghost_get<0>();
@@ -1162,9 +1231,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
// Create a verlet list for each particle
- openfpm::vector<openfpm::vector<size_t>> verlet;
-
- vd.getVerlet(verlet,third_dist);
+ VerletList<3,float,FAST,shift<3,float>> verlet = vd.getVerlet(third_dist);
bool correct = true;
@@ -1179,11 +1246,9 @@ BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
Point<3,float> p = vd.getPos(i);
// for each neighborhood particle
- for (size_t j = 0 ; j < verlet.get(i).size() ; j++)
+ for (size_t j = 0 ; j < verlet.getNNPart(i) ; j++)
{
- auto & NN = verlet.get(i);
-
- Point<3,float> q = vd.getPos(NN.get(j));
+ Point<3,float> q = vd.getPos(verlet.get(i,j));
float dist = p.distance(q);
@@ -1195,7 +1260,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
third_NN++;
}
- correct &= (first_NN == 6);
+ correct &= (first_NN == 7);
correct &= (second_NN == 12);
correct &= (third_NN == 8);
}
@@ -1324,323 +1389,7 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_map_list )
}
}
-///////////////////////// test hilb ///////////////////////////////
-
-BOOST_AUTO_TEST_CASE( vector_dist_reorder_2d_test )
-{
- Vcluster & v_cl = create_vcluster();
-
- // set the seed
- // create the random generator engine
- std::srand(v_cl.getProcessUnitID());
- std::default_random_engine eg;
- std::uniform_real_distribution<float> ud(0.0f, 1.0f);
-
- long int k = 524288 * v_cl.getProcessingUnits();
-
- long int big_step = k / 4;
- big_step = (big_step == 0)?1:big_step;
-
- print_test_v( "Testing 2D vector with hilbert curve reordering k<=",k);
-
- // 2D test
- for ( ; k >= 2 ; k-= decrement(k,big_step) )
- {
- BOOST_TEST_CHECKPOINT( "Testing 2D vector with hilbert curve reordering k=" << k );
-
- Box<2,float> box({0.0,0.0},{1.0,1.0});
-
- // Boundary conditions
- size_t bc[2]={NON_PERIODIC,NON_PERIODIC};
-
- vector_dist<2,float, Point_test<float>, CartDecomposition<2,float> > vd(k,box,bc,Ghost<2,float>(0.0));
-
- auto it = vd.getIterator();
-
- while (it.isNext())
- {
- auto key = it.get();
-
- vd.getPos(key)[0] = ud(eg);
- vd.getPos(key)[1] = ud(eg);
-
- ++it;
- }
-
- vd.map();
-
- // Create first cell list
-
- auto NN1 = vd.getCellList(0.01);
-
- //An order of a curve
- int32_t m = 6;
-
- //Reorder a vector
- vd.reorder(m);
-
- // Create second cell list
- auto NN2 = vd.getCellList(0.01);
-
- //Check equality of cell sizes
- for (size_t i = 0 ; i < NN1.getGrid().size() ; i++)
- {
- size_t n1 = NN1.getNelements(i);
- size_t n2 = NN2.getNelements(i);
-
- BOOST_REQUIRE_EQUAL(n1,n2);
- }
- }
-}
-
-BOOST_AUTO_TEST_CASE( vector_dist_cl_random_vs_hilb_forces_test )
-{
- Vcluster & v_cl = create_vcluster();
-
- if (v_cl.getProcessingUnits() > 32)
- return;
-
- ///////////////////// INPUT DATA //////////////////////
-
- // Dimensionality of the space
- const size_t dim = 3;
- // Cut-off radiuses. Can be put different number of values
- openfpm::vector<float> cl_r_cutoff {0.05};
- // The starting amount of particles (remember that this number is multiplied by number of processors you use for testing)
- size_t cl_k_start = 10000;
- // The lower threshold for number of particles
- size_t cl_k_min = 1000;
- // Ghost part of distributed vector
- double ghost_part = 0.05;
-
- ///////////////////////////////////////////////////////
-
- //For different r_cut
- for (size_t r = 0; r < cl_r_cutoff.size(); r++ )
- {
- //Cut-off radius
- float r_cut = cl_r_cutoff.get(r);
-
- //Number of particles
- size_t k = cl_k_start * v_cl.getProcessingUnits();
-
- std::string str("Testing " + std::to_string(dim) + "D vector's forces (random vs hilb celllist) k<=");
-
- vector_dist_test::print_test_v(str,k);
-
- //For different number of particles
- for (size_t k_int = k ; k_int >= cl_k_min ; k_int/=2 )
- {
- BOOST_TEST_CHECKPOINT( "Testing " << dim << "D vector's forces (random vs hilb celllist) k<=" << k_int );
-
- Box<dim,float> box;
-
- for (size_t i = 0; i < dim; i++)
- {
- box.setLow(i,0.0);
- box.setHigh(i,1.0);
- }
-
- // Boundary conditions
- size_t bc[dim];
-
- for (size_t i = 0; i < dim; i++)
- bc[i] = PERIODIC;
-
- vector_dist<dim,float, aggregate<float[dim]>, CartDecomposition<dim,float> > vd(k_int,box,bc,Ghost<dim,float>(ghost_part));
-
- vector_dist<dim,float, aggregate<float[dim]>, CartDecomposition<dim,float> > vd2(k_int,box,bc,Ghost<dim,float>(ghost_part));
-
- // Initialize dist vectors
- vd_initialize_double<dim>(vd, vd2, v_cl, k_int);
-
- vd.template ghost_get<0>();
- vd2.template ghost_get<0>();
-
- //Get a cell list
-
- auto NN = vd.getCellList(r_cut);
-
- //Calculate forces
-
- calc_forces<dim>(NN,vd,r_cut);
-
- //Get a cell list hilb
-
- auto NN_hilb = vd2.getCellList_hilb(r_cut);
-
- //Calculate forces
- calc_forces_hilb<dim>(NN_hilb,vd2,r_cut);
-
- // Calculate average
- size_t count = 1;
- Point<dim,float> avg;
- for (size_t i = 0 ; i < dim ; i++) {avg.get(i) = 0.0;}
-
- auto it_v2 = vd.getIterator();
- while (it_v2.isNext())
- {
- //key
- vect_dist_key_dx key = it_v2.get();
-
- for (size_t i = 0; i < dim; i++)
- avg.get(i) += fabs(vd.template getProp<0>(key)[i]);
-
- ++count;
- ++it_v2;
- }
-
- for (size_t i = 0 ; i < dim ; i++) {avg.get(i) /= count;}
-
- auto it_v = vd.getIterator();
- while (it_v.isNext())
- {
- //key
- vect_dist_key_dx key = it_v.get();
-
- for (size_t i = 0; i < dim; i++)
- {
- auto a1 = vd.template getProp<0>(key)[i];
- auto a2 = vd2.template getProp<0>(key)[i];
-
- //Check that the forces are (almost) equal
- float per = 0.1;
- if (a1 != 0.0)
- per = fabs(0.1*avg.get(i)/a1);
-
- BOOST_REQUIRE_CLOSE(a1,a2,per);
- }
-
- ++it_v;
- }
- }
- }
-}
-
-BOOST_AUTO_TEST_CASE( vector_dist_cl_random_vs_reorder_forces_test )
-{
- ///////////////////// INPUT DATA //////////////////////
-
- // Dimensionality of the space
- const size_t dim = 3;
- // Cut-off radiuses. Can be put different number of values
- openfpm::vector<float> cl_r_cutoff {0.01};
- // The starting amount of particles (remember that this number is multiplied by number of processors you use for testing)
- size_t cl_k_start = 10000;
- // The lower threshold for number of particles
- size_t cl_k_min = 1000;
- // Ghost part of distributed vector
- double ghost_part = 0.01;
-
- ///////////////////////////////////////////////////////
-
- //For different r_cut
- for (size_t r = 0; r < cl_r_cutoff.size(); r++ )
- {
- Vcluster & v_cl = create_vcluster();
-
- //Cut-off radius
- float r_cut = cl_r_cutoff.get(r);
-
- //Number of particles
- size_t k = cl_k_start * v_cl.getProcessingUnits();
-
- std::string str("Testing " + std::to_string(dim) + "D vector's forces (random vs reorder) k<=");
-
- vector_dist_test::print_test_v(str,k);
-
- //For different number of particles
- for (size_t k_int = k ; k_int >= cl_k_min ; k_int/=2 )
- {
- BOOST_TEST_CHECKPOINT( "Testing " << dim << "D vector's forces (random vs reorder) k<=" << k_int );
-
- Box<dim,float> box;
-
- for (size_t i = 0; i < dim; i++)
- {
- box.setLow(i,0.0);
- box.setHigh(i,1.0);
- }
-
- // Boundary conditions
- size_t bc[dim];
-
- for (size_t i = 0; i < dim; i++)
- bc[i] = PERIODIC;
-
- vector_dist<dim,float, aggregate<float[dim], float[dim]>, CartDecomposition<dim,float> > vd(k_int,box,bc,Ghost<dim,float>(ghost_part));
-
- // Initialize vd
- vd_initialize<dim,decltype(vd)>(vd, v_cl, k_int);
-
- vd.template ghost_get<0>();
-
- //Get a cell list
-
- auto NN1 = vd.getCellList(r_cut);
-
- //Calculate forces '0'
-
- calc_forces<dim>(NN1,vd,r_cut);
-
- //Reorder and get a cell list again
-
- vd.reorder(4);
-
- vd.template ghost_get<0>();
-
- auto NN2 = vd.getCellList(r_cut);
-
- //Calculate forces '1'
- calc_forces<dim,1>(NN2,vd,r_cut);
-
- // Calculate average (For Coverty scan we start from 1)
- size_t count = 1;
- Point<dim,float> avg;
- for (size_t i = 0 ; i < dim ; i++) {avg.get(i) = 0.0;}
-
- auto it_v2 = vd.getIterator();
- while (it_v2.isNext())
- {
- //key
- vect_dist_key_dx key = it_v2.get();
-
- for (size_t i = 0; i < dim; i++)
- avg.get(i) += fabs(vd.template getProp<0>(key)[i]);
-
- ++count;
- ++it_v2;
- }
-
- for (size_t i = 0 ; i < dim ; i++) {avg.get(i) /= count;}
-
- //Test for equality of forces
- auto it_v = vd.getDomainIterator();
-
- while (it_v.isNext())
- {
- //key
- vect_dist_key_dx key = it_v.get();
-
- for (size_t i = 0; i < dim; i++)
- {
- auto a1 = vd.template getProp<0>(key)[i];
- auto a2 = vd.template getProp<1>(key)[i];
-
- //Check that the forces are (almost) equal
- float per = 0.1;
- if (a1 != 0.0)
- per = fabs(0.1*avg.get(i)/a1);
-
- BOOST_REQUIRE_CLOSE(a1,a2,per);
- }
-
- ++it_v;
- }
- }
- }
-}
-
+#include "vector_dist_cell_list_tests.hpp"
BOOST_AUTO_TEST_SUITE_END()