Failed global reductions

example/Vector/3_molecular_dynamic/Makefile

@@ -6,11 +6,13 @@ LDIR =
 
 OBJ = main.o
 OBJ_EXPR = main_expr.o
+OBJ_VL = main_vl.o
+OBJ_VL_SYM = main_vl_sym.o
 
-all: md_dyn md_dyn_expr
+all: md_dyn md_dyn_expr md_dyn_vl md_dyn_vl_sym
 
 %.o: %.cpp
-	$(CC) -O3 -g -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
+	$(CC) -O3 -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
 
 md_dyn: $(OBJ)
 	$(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
@@ -18,11 +20,17 @@ md_dyn: $(OBJ)
 md_dyn_expr: $(OBJ_EXPR)
 	$(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
 
+md_dyn_vl: $(OBJ_VL)
+	$(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
+
+md_dyn_vl_sym: $(OBJ_VL_SYM)
+	$(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
+
 run: all
-	source $$HOME/openfpm_vars; mpirun -np 3 ./md_dyn; mpirun -np 3 ./md_dyn_expr
+	source $$HOME/openfpm_vars; mpirun -np 3 ./md_dyn; mpirun -np 3 ./md_dyn_expr; mpirun -np 3 ./md_dyn_vl; mpirun -np 3 ./md_dyn_vl_sym
 
 .PHONY: clean all run
 
 clean:
-	rm -f *.o *~ core md_dyn md_dyn_expr
+	rm -f *.o *~ core md_dyn md_dyn_expr md_dyn_vl md_dyn_vl_sym
 

example/Vector/3_molecular_dynamic/main.cpp

-
 #include "Vector/vector_dist.hpp"
 #include "Decomposition/CartDecomposition.hpp"
 #include "data_type/aggregate.hpp"
 #include "Plot/GoogleChart.hpp"
 #include "Plot/util.hpp"
+#include "timer.hpp"
 
 /*!
  * \page Vector_3_md Vector 3 molecular dynamic
@@ -48,7 +48,7 @@ constexpr int force = 1;
 
 //! \cond [calc forces] \endcond
 
-void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6)
+void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6, double r_cut2)
 {
 
 //! \cond [calc forces] \endcond
@@ -127,6 +127,9 @@ void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, Ce
 			// take the norm of this vector
 			double rn = norm2(r);
 
+			if (rn > r_cut2)
+			{++Np; continue;};
+
 			// Calculate the force, using pow is slower
 			Point<3,double> f = 24.0*(2.0 *sigma12 / (rn*rn*rn*rn*rn*rn*rn) -  sigma6 / (rn*rn*rn*rn)) * r;
 
@@ -165,8 +168,10 @@ void calc_forces(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, Ce
 
 //! \cond [calc energy] \endcond
 
-double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6)
+double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd, CellList<3, double, FAST, shift<3, double> > & NN, double sigma12, double sigma6, double r_cut2)
 {
+        double rc = r_cut2;
+        double shift = 2.0 * ( sigma12 / (rc*rc*rc*rc*rc*rc) - sigma6 / ( rc*rc*rc) );
 
 //! \cond [calc energy] \endcond
 
@@ -185,7 +190,7 @@ double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd,
 	//! \cond [up cell ene] \endcond
 
 	double E = 0.0;
-	vd.updateCellList(NN);
+//	vd.updateCellList(NN);
 
 	//! \cond [up cell ene] \endcond
 
@@ -236,8 +241,11 @@ double calc_energy(vector_dist<3,double, aggregate<double[3],double[3]> > & vd,
 			// take the normalized direction
 			double rn = norm2(xp - xq);
 
+			if (rn > r_cut2)
+			{++Np;continue;}
+
 			// potential energy (using pow is slower)
-			E += 2.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) ) - shift;
 
 			// Next neighborhood
 			++Np;
@@ -414,13 +422,16 @@ int main(int argc, char* argv[])
 	 *
 	 */
 
+	timer tsim;
+	tsim.start();
+
 	//! \cond [md steps] \endcond
 
 	// Get the Cell list structure
 	auto NN = vd.getCellList(r_cut);
 
 	// calculate forces
-	calc_forces(vd,NN,sigma12,sigma6);
+	calc_forces(vd,NN,sigma12,sigma6,r_cut*r_cut);
 	unsigned long int f = 0;
 
 	// MD time stepping
@@ -452,7 +463,7 @@ int main(int argc, char* argv[])
 		vd.template ghost_get<>();
 
 		// calculate forces or a(tn + 1) Step 2
-		calc_forces(vd,NN,sigma12,sigma6);
+		calc_forces(vd,NN,sigma12,sigma6,r_cut*r_cut);
 
 
 		// Integrate the velocity Step 3
@@ -481,7 +492,7 @@ int main(int argc, char* argv[])
 			vd.ghost_get<>();
 
 			// We calculate the energy
-			double energy = calc_energy(vd,NN,sigma12,sigma6);
+			double energy = calc_energy(vd,NN,sigma12,sigma6,r_cut*r_cut);
 			auto & vcl = create_vcluster();
 			vcl.sum(energy);
 			vcl.execute();
@@ -501,6 +512,9 @@ int main(int argc, char* argv[])
 
 	//! \cond [md steps] \endcond
 
+	tsim.stop();
+	std::cout << "Time: " << tsim.getwct() << std::endl;
+
 	/*!
 	 * \page Vector_3_md Vector 3 molecular dynamic
 	 *

example/Vector/3_molecular_dynamic/main_expr.cpp

 #include "Vector/vector_dist.hpp"
 #include "Plot/GoogleChart.hpp"
 #include "Operators/Vector/vector_dist_operators.hpp"
+#include "timer.hpp"
 
 constexpr int velocity = 0;
 constexpr int force = 1;
 
 struct ln_potential
 {
-	double sigma12,sigma6;
+	double sigma12,sigma6,r_cut2,shift;
 
-	ln_potential(double sigma12_, double sigma6_) {sigma12 = sigma12_, sigma6 = sigma6_;}
+	ln_potential(double sigma12_, double sigma6_, double r_cut2_, double shift_) {sigma12 = sigma12_; sigma6 = sigma6_; r_cut2 = r_cut2_;shift = shift_;}
 
 	Point<2,double> value(const Point<3,double> & xp, const Point<3,double> xq)
 	{
 		double rn = norm2(xp - xq);
+		if (rn >= r_cut2)	return 0.0;
 
-		Point<2,double> E({2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ),
-						 0.0});
+		Point<2,double> E({2.0 * ( sigma12 / (rn*rn*rn*rn*rn*rn) - sigma6 / ( rn*rn*rn) ) - shift,0.0});
 
 		return E;
 	}
@@ -24,24 +25,28 @@ struct ln_potential
 
 struct ln_force
 {
-	double sigma12,sigma6;
+	double sigma12,sigma6,r_cut2;
 
-	ln_force(double sigma12_, double sigma6_) {sigma12 = sigma12_; sigma6 = sigma6_;}
+	ln_force(double sigma12_, double sigma6_, double r_cut2_) {sigma12 = sigma12_; sigma6 = sigma6_;r_cut2 = r_cut2_;}
 
 	Point<3,double> value(const Point<3,double> & xp, const Point<3,double> xq)
 	{
 		Point<3,double> r = xp - xq;
 		double rn = norm2(r);
 
+		if (rn > r_cut2)	return 0.0;
+
 		return 24.0*(2.0 * sigma12 / (rn*rn*rn*rn*rn*rn*rn) -  sigma6 / (rn*rn*rn*rn)) * r;
 	}
 };
 
 int main(int argc, char* argv[])
 {
-	double dt = 0.0005, sigma = 0.1, r_cut = 0.3;
+	double dt = 0.0005, sigma = 0.1, r_cut = 3.0*sigma;
 
 	double sigma6 = pow(sigma,6), sigma12 = pow(sigma,12);
+	double rc2 = r_cut * r_cut;
+	double shift = 2.0 * ( sigma12 / (rc2*rc2*rc2*rc2*rc2*rc2) - sigma6 / ( rc2*rc2*rc2) );
 
 	openfpm::vector<double> x;
 	openfpm::vector<openfpm::vector<double>> y;
@@ -55,8 +60,8 @@ int main(int argc, char* argv[])
 	size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
 	Ghost<3,float> ghost(r_cut);
 
-	ln_force lf(sigma12,sigma6);
-	ln_potential lp(sigma12,sigma6);
+	ln_force lf(sigma12,sigma6,r_cut*r_cut);
+	ln_potential lp(sigma12,sigma6,r_cut*r_cut,shift);
 
 	vector_dist<3,double, aggregate<Point<3,double>,Point<3,double>> > vd(0,box,bc,ghost);
 
@@ -82,6 +87,9 @@ int main(int argc, char* argv[])
 	v_force = 0;
 	v_velocity = 0;
 
+	timer tsim;
+	tsim.start();
+
 	auto NN = vd.getCellList(r_cut);
 
 	vd.updateCellList(NN);
@@ -109,6 +117,7 @@ int main(int argc, char* argv[])
 			vd.write("particles_",f);
 			vd.ghost_get<>();
 
+			vd.updateCellList(NN);
 			Point<2,double> E = rsum(applyKernel_in_sim(vd,NN,lp) + (v_velocity * v_velocity)/2.0,vd).get();
 
 			vcl.sum(E.get(0));vcl.sum(E.get(1));
@@ -125,6 +134,9 @@ int main(int argc, char* argv[])
 		}
 	}
 
+	tsim.stop();
+	std::cout << "Time: " << tsim.getwct() << std::endl;
+
 	GCoptions options;
 	options.title = std::string("Energy with time");
 	options.yAxis = std::string("Energy");

openfpm_data @ 9301fe45

-Subproject commit f075dd5da60f5170e01a3265a88881e98fc3b7e4
+Subproject commit 9301fe459ca2ef81115c642ee7a54c983ddfd6f1

src/Makefile.am

@@ -9,7 +9,7 @@ nobase_include_HEADERS = Decomposition/CartDecomposition.hpp Decomposition/CartD
          Decomposition/nn_processor.hpp Decomposition/ie_loc_ghost.hpp Decomposition/ORB.hpp \
          Graph/CartesianGraphFactory.hpp \
          Grid/grid_dist_id.hpp Grid/grid_dist_id_iterator_dec.hpp Grid/grid_dist_util.hpp  Grid/grid_dist_id_iterator_sub.hpp Grid/grid_dist_id_iterator.hpp Grid/grid_dist_key.hpp Grid/staggered_dist_grid.hpp Grid/staggered_dist_grid_util.hpp Grid/staggered_dist_grid_copy.hpp \
-         Vector/vector_dist.hpp Vector/vector_dist_ofb.hpp Vector/vector_dist_iterator.hpp Vector/vector_dist_key.hpp \
+         Vector/vector_dist_comm.hpp Vector/vector_dist.hpp Vector/vector_dist_ofb.hpp Vector/vector_dist_iterator.hpp Vector/vector_dist_key.hpp \
          config/config.h \
          example.mk \
          Decomposition/Distribution/metis_util.hpp Decomposition/Distribution/parmetis_dist_util.hpp  Decomposition/Distribution/parmetis_util.hpp Decomposition/Distribution/MetisDistribution.hpp Decomposition/Distribution/ParMetisDistribution.hpp Decomposition/Distribution/DistParMetisDistribution.hpp  dec_optimizer.hpp SubdomainGraphNodes.hpp \

src/Vector/vector_dist_unit_test.hpp

@@ -1389,7 +1389,192 @@ BOOST_AUTO_TEST_CASE( vector_dist_periodic_map_list )
 	}
 }
 
+
+BOOST_AUTO_TEST_CASE( vector_dist_ghost_with_ghost_buffering )
+{
+	Vcluster & v_cl = create_vcluster();
+
+	if (v_cl.getProcessingUnits() > 3)
+		return;
+
+    // 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("Testing 3D periodic vector with ghost buffering k=",k);
+	BOOST_TEST_CHECKPOINT( "Testing 3D periodic with ghost buffering k=" << k );
+
+	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
+	Ghost<3,float> ghost(0.1);
+
+	typedef  aggregate<float> part_prop;
+
+	// Distributed vector
+	vector_dist<3,float, part_prop > vd(k,box,bc,ghost);
+
+	auto it = vd.getIterator();
+
+	while (it.isNext())
+	{
+		auto key = it.get();
+
+		vd.getPos(key)[0] = ud(eg);
+		vd.getPos(key)[1] = ud(eg);
+		vd.getPos(key)[2] = ud(eg);
+
+		// Fill some properties randomly
+
+		vd.getProp<0>(key) = 0.0;
+
+		++it;
+	}
+
+	vd.map();
+
+	// sync the ghost
+	vd.ghost_get<0>();
+
+	openfpm::vector<size_t> list_idx;
+	openfpm::vector<size_t> list_idx2;
+
+	auto it3 = vd.getGhostIterator();
+	while (it3.isNext())
+	{
+		auto key = it3.get();
+
+		list_idx.add(key.getKey());
+
+		++it3;
+	}
+
+	list_idx.sort();
+
+	for (size_t i = 0 ; i < 10 ; i++)
+	{
+		auto it = vd.getDomainIterator();
+
+		while (it.isNext())
+		{
+			auto key = it.get();
+
+			vd.getPos(key)[0] = ud(eg);
+			vd.getPos(key)[1] = ud(eg);
+			vd.getPos(key)[2] = ud(eg);
+
+			// Fill some properties randomly
+
+			vd.getProp<0>(key) = i;
+
+			++it;
+		}
+
+		vd.ghost_get<0>(SKIP_LABELLING);
+
+		list_idx2.clear();
+		auto it2 = vd.getGhostIterator();
+		bool ret = true;
+
+		while (it2.isNext())
+		{
+			auto key = it2.get();
+
+			list_idx2.add(key.getKey());
+			ret &= vd.getProp<0>(key) == i;
+
+			++it2;
+		}
+
+		BOOST_REQUIRE_EQUAL(ret,true);
+		BOOST_REQUIRE_EQUAL(list_idx.size(),list_idx2.size());
+
+		list_idx2.sort();
+
+		ret = true;
+		for (size_t i = 0 ; i < list_idx.size() ; i++)
+			ret &= list_idx.get(i) == list_idx2.get(i);
+
+		BOOST_REQUIRE_EQUAL(ret,true);
+	}
+
+	vd.map();
+	vd.ghost_get<0>();
+
+	list_idx.clear();
+
+	auto it4 = vd.getGhostIterator();
+	while (it4.isNext())
+	{
+		auto key = it4.get();
+
+		list_idx.add(key.getKey());
+
+		++it4;
+	}
+
+	list_idx.sort();
+
+	for (size_t i = 0 ; i < 10 ; i++)
+	{
+		auto it = vd.getDomainIterator();
+
+		while (it.isNext())
+		{
+			auto key = it.get();
+
+			vd.getPos(key)[0] = ud(eg);
+			vd.getPos(key)[1] = ud(eg);
+			vd.getPos(key)[2] = ud(eg);
+
+			// Fill some properties randomly
+
+			vd.getProp<0>(key) = i;
+
+			++it;
+		}
+
+		vd.ghost_get<0>(SKIP_LABELLING);
+
+		list_idx2.clear();
+		auto it2 = vd.getGhostIterator();
+		bool ret = true;
+
+		while (it2.isNext())
+		{
+			auto key = it2.get();
+
+			list_idx2.add(key.getKey());
+			ret &= vd.getProp<0>(key) == i;
+
+			++it2;
+		}
+
+		BOOST_REQUIRE_EQUAL(ret,true);
+		BOOST_REQUIRE_EQUAL(list_idx.size(),list_idx2.size());
+
+		list_idx2.sort();
+
+		ret = true;
+		for (size_t i = 0 ; i < list_idx.size() ; i++)
+			ret &= list_idx.get(i) == list_idx2.get(i);
+
+		BOOST_REQUIRE_EQUAL(ret,true);
+	}
+}
+
 #include "vector_dist_cell_list_tests.hpp"
+#include "vector_dist_NN_tests.hpp"
 
 BOOST_AUTO_TEST_SUITE_END()
 

src/pdata_performance.hpp

+/*
+ * performance.hpp
+ *
+ *  Created on: Mar 9, 2016
+ *      Author: yaroslav
+ */
+
+#ifndef SRC_PDATA_PERFORMANCE_HPP_
+#define SRC_PDATA_PERFORMANCE_HPP_
+
+BOOST_AUTO_TEST_SUITE( performance )
+
+#include "Vector/vector_dist_verlet_performance_tests.hpp"
+#include "Vector/vector_dist_cl_performance_tests.hpp"
+#include "Vector/vector_dist_cl_hilb_performance_tests.hpp"
+
+
+BOOST_AUTO_TEST_SUITE_END()
+
+#endif /* SRC_PDATA_PERFORMANCE_HPP_ */