vector_dist_unit_test.hpp 39.5 KB
Newer Older
incardon's avatar
incardon committed
1 2 3 4 5 6 7 8 9 10
/*
 * vector_dist_unit_test.hpp
 *
 *  Created on: Mar 6, 2015
 *      Author: Pietro Incardona
 */

#ifndef VECTOR_DIST_UNIT_TEST_HPP_
#define VECTOR_DIST_UNIT_TEST_HPP_

11 12
#include "config.h"

incardon's avatar
incardon committed
13
#include <random>
incardon's avatar
incardon committed
14
#include "Vector/vector_dist.hpp"
15
#include "data_type/aggregate.hpp"
incardon's avatar
incardon committed
16
#include "Vector/performance/vector_dist_performance_common.hpp"
incardon's avatar
incardon committed
17

18 19 20 21 22 23
/*! \brief Count the total number of particles
 *
 * \param vd distributed vector
 * \param bc boundary conditions
 *
 */
incardon's avatar
incardon committed
24
template<unsigned int dim, template <typename> class layout> size_t total_n_part_lc(vector_dist<dim,float, Point_test<float>,typename layout<Point_test<float>>::type, layout, CartDecomposition<dim,float> > & vd, size_t (& bc)[dim])
25 26 27 28 29
{
	Vcluster & v_cl = vd.getVC();
	auto it2 = vd.getDomainIterator();
	const CartDecomposition<3,float> & ct = vd.getDecomposition();

30 31
	bool noOut = true;

32 33 34 35 36
	size_t cnt = 0;
	while (it2.isNext())
	{
		auto key = it2.get();

Pietro Incardona's avatar
Pietro Incardona committed
37
		noOut &= ct.isLocal(vd.getPos(key));
38 39 40 41 42 43

		cnt++;

		++it2;
	}

44 45
	BOOST_REQUIRE_EQUAL(noOut,true);

46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
	//
	v_cl.sum(cnt);
	v_cl.execute();

	return cnt;
}

/*! \brief Count local and non local
 *
 * \param vd distributed vector
 * \param it iterator
 * \param bc boundary conditions
 * \param box domain box
 * \param dom_ext domain + ghost box
 * \param l_cnt local particles counter
 * \param nl_cnt non local particles counter
 * \param n_out out of domain + ghost particles counter
 *
 */
65
template<unsigned int dim,typename vector_dist> inline void count_local_n_local(vector_dist & vd, vector_dist_iterator & it, size_t (& bc)[dim] , Box<dim,float> & box, Box<dim,float> & dom_ext, size_t & l_cnt, size_t & nl_cnt, size_t & n_out)
66 67 68 69 70 71 72
{
	const CartDecomposition<dim,float> & ct = vd.getDecomposition();

	while (it.isNext())
	{
		auto key = it.get();
		// Check if it is in the domain
Pietro Incardona's avatar
Pietro Incardona committed
73
		if (box.isInsideNP(vd.getPos(key)) == true)
74 75
		{
			// Check if local
Pietro Incardona's avatar
Pietro Incardona committed
76
			if (ct.isLocalBC(vd.getPos(key),bc) == true)
77 78 79 80 81 82 83 84 85
				l_cnt++;
			else
				nl_cnt++;
		}
		else
		{
			nl_cnt++;
		}

incardon's avatar
incardon committed
86 87
		Point<dim,float> xp = vd.getPos(key);

88
		// Check that all particles are inside the Domain + Ghost part
incardon's avatar
incardon committed
89
		if (dom_ext.isInside(xp) == false)
90 91 92 93 94 95
				n_out++;

		++it;
	}
}

incardon's avatar
incardon committed
96 97
BOOST_AUTO_TEST_SUITE( vector_dist_test )

98 99
void print_test(std::string test, size_t sz)
{
100
	if (create_vcluster().getProcessUnitID() == 0)
101 102 103
		std::cout << test << " " << sz << "\n";
}

104
void Test2D_ghost(Box<2,float> & box)
105 106
{
	// Communication object
107
	Vcluster & v_cl = create_vcluster();
108 109 110

	typedef Point_test<float> p;

111
	// Get the default minimum number of sub-sub-domain per processor (granularity of the decomposition)
incardon's avatar
incardon committed
112
	size_t n_sub = 64 * v_cl.getProcessingUnits();
113 114
	// Convert the request of having a minimum n_sub number of sub-sub domain into grid decompsition of the space
	size_t sz = CartDecomposition<2,float>::getDefaultGrid(n_sub);
incardon's avatar
incardon committed
115

116 117
	//! [Create a vector of elements distributed on a grid like way]

118
	size_t g_div[]= {sz,sz};
119

120 121 122
	// number of particles
	size_t np = sz * sz;

123
	// Calculate the number of elements this processor is going to obtain
124 125 126 127 128 129 130 131
	size_t p_np = np / v_cl.getProcessingUnits();

	// Get non divisible part
	size_t r = np % v_cl.getProcessingUnits();

	// Get the offset
	size_t offset = v_cl.getProcessUnitID() * p_np + std::min(v_cl.getProcessUnitID(),r);

132
	// Distribute the remain elements
133 134
	if (v_cl.getProcessUnitID() < r)
		p_np++;
135 136 137 138 139

	// Create a grid info
	grid_sm<2,void> g_info(g_div);

	// Calculate the grid spacing
140
	Point<2,float> spacing = box.getP2() - box.getP1();
141 142 143
	spacing = spacing / g_div;

	// middle spacing
Pietro Incardona's avatar
Pietro Incardona committed
144
	Point<2,float> m_spacing = spacing / 2.0;
145

incardon's avatar
incardon committed
146 147 148
	// set the ghost based on the radius cut off (make just a little bit smaller than the spacing)
	Ghost<2,float> g(spacing.get(0) - spacing .get(0) * 0.0001);

149 150 151
	// Boundary conditions
	size_t bc[2]={NON_PERIODIC,NON_PERIODIC};

152
	// Vector of particles
153
	vector_dist<2,float, Point_test<float> > vd(g_info.size(),box,bc,g);
154 155 156

	// size_t
	size_t cobj = 0;
157

158 159
	grid_key_dx_iterator_sp<2> it(g_info,offset,offset+p_np-1);
	auto v_it = vd.getIterator();
160

161
	while (v_it.isNext() && it.isNext())
162
	{
163 164
		auto key = it.get();
		auto key_v = v_it.get();
165 166 167

		// set the particle position

168 169
		vd.getPos(key_v)[0] = key.get(0) * spacing[0] + m_spacing[0] + box.getLow(0);
		vd.getPos(key_v)[1] = key.get(1) * spacing[1] + m_spacing[1] + box.getLow(1);
170

171 172 173
		cobj++;

		++v_it;
174 175 176
		++it;
	}

177 178 179
	//! [Create a vector of elements distributed on a grid like way]

	// Both iterators must signal the end, and the number of elements in the vector, must the equal to the
180 181 182 183 184
	// predicted one
	BOOST_REQUIRE_EQUAL(v_it.isNext(),false);
	BOOST_REQUIRE_EQUAL(it.isNext(),false);
	BOOST_REQUIRE_EQUAL(cobj,p_np);

185 186
	//! [Redistribute the particles and sync the ghost properties]

187 188 189
	// redistribute the particles according to the decomposition
	vd.map();

190
	auto v_it2 = vd.getIterator();
191

192
	while (v_it2.isNext())
193
	{
194
		auto key = v_it2.get();
195 196

		// fill with the processor ID where these particle live
incardon's avatar
incardon committed
197
		vd.getProp<p::s>(key) = vd.getPos(key)[0] + vd.getPos(key)[1] * 16.0f;
198 199 200
		vd.getProp<p::v>(key)[0] = v_cl.getProcessUnitID();
		vd.getProp<p::v>(key)[1] = v_cl.getProcessUnitID();
		vd.getProp<p::v>(key)[2] = v_cl.getProcessUnitID();
201

202
		++v_it2;
incardon's avatar
incardon committed
203
	}
incardon's avatar
incardon committed
204 205

	// do a ghost get
206
	vd.ghost_get<p::s,p::v>();
207

208
	//! [Redistribute the particles and sync the ghost properties]
incardon's avatar
incardon committed
209

210 211 212 213
	// Get the decomposition
	const auto & dec = vd.getDecomposition();

	// Get the ghost external boxes
214
	openfpm::vector<size_t> vb(dec.getNEGhostBox());
215 216 217 218

	// Get the ghost iterator
	auto g_it = vd.getGhostIterator();

219 220
	size_t n_part = 0;

221 222 223 224 225 226
	// Check if the ghost particles contain the correct information
	while (g_it.isNext())
	{
		auto key = g_it.get();

		// Check the received data
incardon's avatar
incardon committed
227
		BOOST_REQUIRE_EQUAL(vd.getPos(key)[0] + vd.getPos(key)[1] * 16.0f,vd.getProp<p::s>(key));
228 229 230

		bool is_in = false;
		size_t b = 0;
231
		size_t lb = 0;
232

233
		// check if the received data are in one of the ghost boxes
234
		for ( ; b < dec.getNEGhostBox() ; b++)
235
		{
incardon's avatar
incardon committed
236 237 238
			Point<2,float> xp = vd.getPos(key);

			if (dec.getEGhostBox(b).isInside(xp) == true )
239 240 241 242 243 244 245
			{
				is_in = true;

				// Add
				vb.get(b)++;
				lb = b;
			}
246 247 248 249
		}
		BOOST_REQUIRE_EQUAL(is_in,true);

		// Check that the particle come from the correct processor
250
		BOOST_REQUIRE_EQUAL(vd.getProp<p::v>(key)[0],dec.getEGhostBoxProcessor(lb));
251

252
		n_part++;
253 254 255
		++g_it;
	}

256 257 258 259
	if (v_cl.getProcessingUnits() > 1)
	{
		BOOST_REQUIRE(n_part != 0);
	}
260

261
    CellDecomposer_sm<2,float,shift<2,float>> cd(SpaceBox<2,float>(box),g_div,0);
262 263 264 265

	for (size_t i = 0 ; i < vb.size() ; i++)
	{
		// Calculate how many particle should be in the box
266
		size_t n_point = cd.getGridPoints(dec.getEGhostBox(i)).getVolumeKey();
267

268
		BOOST_REQUIRE_EQUAL(n_point,vb.get(i));
269
	}
incardon's avatar
incardon committed
270 271
}

272 273 274 275 276 277 278 279 280
BOOST_AUTO_TEST_CASE( vector_dist_ghost )
{
	Box<2,float> box({0.0,0.0},{1.0,1.0});
	Test2D_ghost(box);

	Box<2,float> box2({-1.0,-1.0},{2.5,2.5});
	Test2D_ghost(box2);
}

281 282
void print_test_v(std::string test, size_t sz)
{
283
	if (create_vcluster().getProcessUnitID() == 0)
284 285 286
		std::cout << test << " " << sz << "\n";
}

incardon's avatar
incardon committed
287 288 289 290 291 292 293 294 295 296 297 298 299 300
long int decrement(long int k, long int step)
{
	if (k <= 32)
	{
		return 1;
	}
	else if (k - 2*step+1 <= 0)
	{
		return k - 32;
	}
	else
		return step;
}

301
BOOST_AUTO_TEST_CASE( vector_dist_iterator_test_use_2d )
incardon's avatar
incardon committed
302
{
303
	Vcluster & v_cl = create_vcluster();
incardon's avatar
incardon committed
304 305 306 307 308 309 310

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

311 312 313
#ifdef TEST_COVERAGE_MODE
    long int k = 24288 * v_cl.getProcessingUnits();
#else
314
    long int k = 524288 * v_cl.getProcessingUnits();
315
#endif
incardon's avatar
incardon committed
316

317
	long int big_step = k / 4;
318
	big_step = (big_step == 0)?1:big_step;
incardon's avatar
incardon committed
319

320 321
	print_test_v( "Testing 2D vector k<=",k);

322
	// 2D test
incardon's avatar
incardon committed
323
	for ( ; k >= 2 ; k-= decrement(k,big_step) )
incardon's avatar
incardon committed
324
	{
325
		BOOST_TEST_CHECKPOINT( "Testing 2D vector k=" << k );
326 327 328

		//! [Create a vector of random elements on each processor 2D]

329
		Box<2,float> box({0.0,0.0},{1.0,1.0});
330 331 332 333

		// Boundary conditions
		size_t bc[2]={NON_PERIODIC,NON_PERIODIC};

334
		vector_dist<2,float, Point_test<float> > vd(k,box,bc,Ghost<2,float>(0.0));
incardon's avatar
incardon committed
335

336
		auto it = vd.getIterator();
incardon's avatar
incardon committed
337

338 339 340 341
		while (it.isNext())
		{
			auto key = it.get();

Pietro Incardona's avatar
Pietro Incardona committed
342 343
			vd.getPos(key)[0] = ud(eg);
			vd.getPos(key)[1] = ud(eg);
344 345 346 347 348 349

			++it;
		}

		vd.map();

350 351
		//! [Create a vector of random elements on each processor 2D]

352 353 354
		// Check if we have all the local particles
		size_t cnt = 0;
		const CartDecomposition<2,float> & ct = vd.getDecomposition();
355
		auto it2 = vd.getIterator();
356

357
		while (it2.isNext())
358
		{
359
			auto key = it2.get();
360 361

			// Check if local
Pietro Incardona's avatar
Pietro Incardona committed
362
			BOOST_REQUIRE_EQUAL(ct.isLocal(vd.getPos(key)),true);
363 364 365

			cnt++;

366
			++it2;
367 368 369 370 371
		}

		//
		v_cl.sum(cnt);
		v_cl.execute();
372
		BOOST_REQUIRE_EQUAL((long int)cnt,k);
incardon's avatar
incardon committed
373
	}
374
}
incardon's avatar
incardon committed
375

376 377
BOOST_AUTO_TEST_CASE( vector_dist_iterator_test_use_3d )
{
378
	Vcluster & v_cl = create_vcluster();
incardon's avatar
incardon committed
379

380 381 382 383 384 385
    // 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);

386 387 388
#ifdef TEST_COVERAGE_MODE
    long int k = 24288 * v_cl.getProcessingUnits();
#else
389
    long int k = 524288 * v_cl.getProcessingUnits();
390
#endif
391

392
	long int big_step = k / 4;
393 394
	big_step = (big_step == 0)?1:big_step;

395 396
	print_test_v( "Testing 3D vector k<=",k);

397
	// 3D test
incardon's avatar
incardon committed
398
	for ( ; k >= 2 ; k-= decrement(k,big_step) )
incardon's avatar
incardon committed
399
	{
400
		BOOST_TEST_CHECKPOINT( "Testing 3D vector k=" << k );
401 402 403

		//! [Create a vector of random elements on each processor 3D]

404
		Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
405 406 407 408

		// Boundary conditions
		size_t bc[3]={NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};

409
		vector_dist<3,float, Point_test<float> > vd(k,box,bc,Ghost<3,float>(0.0));
incardon's avatar
incardon committed
410

411
		auto it = vd.getIterator();
incardon's avatar
incardon committed
412

413 414 415
		while (it.isNext())
		{
			auto key = it.get();
incardon's avatar
incardon committed
416

Pietro Incardona's avatar
Pietro Incardona committed
417 418 419
			vd.getPos(key)[0] = ud(eg);
			vd.getPos(key)[1] = ud(eg);
			vd.getPos(key)[2] = ud(eg);
420 421 422 423 424 425

			++it;
		}

		vd.map();

426 427
		//! [Create a vector of random elements on each processor 3D]

428 429 430
		// Check if we have all the local particles
		size_t cnt = 0;
		const CartDecomposition<3,float> & ct = vd.getDecomposition();
431
		auto it2 = vd.getIterator();
432

433
		while (it2.isNext())
434
		{
435
			auto key = it2.get();
436 437

			// Check if local
Pietro Incardona's avatar
Pietro Incardona committed
438
			BOOST_REQUIRE_EQUAL(ct.isLocal(vd.getPos(key)),true);
439 440 441

			cnt++;

442
			++it2;
443 444 445 446 447
		}

		//
		v_cl.sum(cnt);
		v_cl.execute();
448
		BOOST_REQUIRE_EQUAL(cnt,(size_t)k);
449
	}
incardon's avatar
incardon committed
450 451
}

Pietro Incardona's avatar
Pietro Incardona committed
452 453 454 455 456 457 458 459 460 461 462

BOOST_AUTO_TEST_CASE( vector_dist_iterator_fixed_dec_3d )
{
	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);

463 464 465
#ifdef TEST_COVERAGE_MODE
    long int k = 2428 * v_cl.getProcessingUnits();
#else
Pietro Incardona's avatar
Pietro Incardona committed
466
    long int k = 52428 * v_cl.getProcessingUnits();
467
#endif
Pietro Incardona's avatar
Pietro Incardona committed
468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529

	long int big_step = k / 4;
	big_step = (big_step == 0)?1:big_step;

	print_test_v( "Testing 3D vector copy decomposition k<=",k);

	// 3D test
	for ( ; k >= 2 ; k-= decrement(k,big_step) )
	{
		BOOST_TEST_CHECKPOINT( "Testing 3D vector copy decomposition k=" << k );

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

		vector_dist<3,float, aggregate<double,double> > vd(k,box,bc,Ghost<3,float>(0.05));
		vector_dist<3,float, aggregate<double,double> > vd2(vd.getDecomposition(),k);

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

			vd2.getPos(key)[0] = vd.getPos(key)[0];
			vd2.getPos(key)[1] = vd.getPos(key)[1];
			vd2.getPos(key)[2] = vd.getPos(key)[2];

			++it;
		}

		vd.map();
		vd2.map();

		vd.ghost_get();
		vd2.ghost_get();

		auto NN = vd.getCellList(0.05);
		auto NN2 = vd2.getCellList(0.05);

		cross_calc<3,0>(NN,NN2,vd,vd2);
		cross_calc<3,1>(NN,NN,vd,vd);


		auto it3 = vd.getIterator();

		while (it3.isNext())
		{
			auto key = it3.get();

			BOOST_REQUIRE_EQUAL(vd.getProp<0>(key),vd.getProp<1>(key));

			++it3;
		}
	}
}

530 531
BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_use_2d )
{
532
	Vcluster & v_cl = create_vcluster();
533 534 535 536 537 538 539

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

540 541 542
#ifdef TEST_COVERAGE_MODE
    long int k = 24288 * v_cl.getProcessingUnits();
#else
543
    long int k = 524288 * v_cl.getProcessingUnits();
544
#endif
545

546
	long int big_step = k / 4;
547 548 549 550 551 552 553 554 555 556 557 558 559 560
	big_step = (big_step == 0)?1:big_step;

	print_test_v( "Testing 2D periodic vector k<=",k);

	// 2D test
	for ( ; k >= 2 ; k-= decrement(k,big_step) )
	{
		BOOST_TEST_CHECKPOINT( "Testing 2D periodic vector k=" << k );

		Box<2,float> box({0.0,0.0},{1.0,1.0});

		// Boundary conditions
		size_t bc[2]={PERIODIC,PERIODIC};

561
		// factor
562
		float factor = pow(create_vcluster().getProcessingUnits()/2.0f,1.0f/3.0f);
563

564
		// ghost
565
		Ghost<2,float> ghost(0.01 / factor);
566

567
		// ghost2 (a little bigger because of round off error)
568
		Ghost<2,float> ghost2(0.05001 / factor);
569 570

		// Distributed vector
571
		vector_dist<2,float, Point_test<float> > vd(k,box,bc,ghost);
572 573 574 575 576 577 578

		auto it = vd.getIterator();

		while (it.isNext())
		{
			auto key = it.get();

Pietro Incardona's avatar
Pietro Incardona committed
579 580
			vd.getPos(key)[0] = ud(eg);
			vd.getPos(key)[1] = ud(eg);
581 582 583 584 585 586

			++it;
		}

		vd.map();

587
		// sync the ghost, only the property zero
588 589
		vd.ghost_get<0>();

590
		// Domain + ghost box
591
		Box<2,float> dom_ext = box;
592
		dom_ext.enlarge(ghost2);
593

594 595 596 597
		// Iterate on all particles domain + ghost
		size_t l_cnt = 0;
		size_t nl_cnt = 0;
		size_t n_out = 0;
598 599


600
		auto it2 = vd.getIterator();
601
		count_local_n_local<2,vector_dist<2,float, Point_test<float> >>(vd,it2,bc,box,dom_ext,l_cnt,nl_cnt,n_out);
602

603 604
		// No particles should be out of domain + ghost
		BOOST_REQUIRE_EQUAL(n_out,0ul);
605

606
		// Ghost must populated because we synchronized them
607
		if (k > 524288)
608
		{
609
			BOOST_REQUIRE(nl_cnt != 0);
610 611
			BOOST_REQUIRE(l_cnt > nl_cnt);
		}
612

613
		// Sum all the particles inside the domain
614 615
		v_cl.sum(l_cnt);
		v_cl.execute();
616 617

		// count that they are equal to the initial total number
618
		BOOST_REQUIRE_EQUAL((long int)l_cnt,k);
619 620 621 622 623 624

		l_cnt = 0;
		nl_cnt = 0;

		// Iterate only on the ghost particles
		auto itg = vd.getGhostIterator();
625
		count_local_n_local<2,vector_dist<2,float, Point_test<float> > >(vd,itg,bc,box,dom_ext,l_cnt,nl_cnt,n_out);
626 627 628 629 630 631

		// No particle on the ghost must be inside the domain
		BOOST_REQUIRE_EQUAL(l_cnt,0ul);

		// Ghost must be populated
		if (k > 524288)
632
		{
633
			BOOST_REQUIRE(nl_cnt != 0);
634
		}
635 636 637 638 639
	}
}

BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_use_3d )
{
640
	Vcluster & v_cl = create_vcluster();
641 642 643 644 645 646 647

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

648 649 650
#ifdef TEST_COVERAGE_MODE
    long int k = 24288 * v_cl.getProcessingUnits();
#else
651
    long int k = 524288 * v_cl.getProcessingUnits();
652
#endif
653

654
	long int big_step = k / 4;
655 656 657 658 659 660 661 662 663 664 665 666 667 668
	big_step = (big_step == 0)?1:big_step;

	print_test_v( "Testing 3D periodic vector k<=",k);

	// 3D test
	for ( ; k >= 2 ; k-= decrement(k,big_step) )
	{
		BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector 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};

669
		// factor
670
		float factor = pow(create_vcluster().getProcessingUnits()/2.0f,1.0f/3.0f);
671

672
		// ghost
673
		Ghost<3,float> ghost(0.05 / factor);
674

675
		// ghost2 (a little bigger because of round off error)
676
		Ghost<3,float> ghost2(0.05001 / factor);
677 678

		// Distributed vector
Pietro Incardona's avatar
Pietro Incardona committed
679
		vector_dist<3,float, Point_test<float> > vd(k,box,bc,ghost);
680 681 682 683 684 685 686

		auto it = vd.getIterator();

		while (it.isNext())
		{
			auto key = it.get();

Pietro Incardona's avatar
Pietro Incardona committed
687 688 689
			vd.getPos(key)[0] = ud(eg);
			vd.getPos(key)[1] = ud(eg);
			vd.getPos(key)[2] = ud(eg);
690 691 692 693 694 695 696 697 698

			++it;
		}

		vd.map();

		// sync the ghost
		vd.ghost_get<0>();

699 700 701 702 703
		// Domain + ghost
		Box<3,float> dom_ext = box;
		dom_ext.enlarge(ghost2);

		// Iterate on all particles domain + ghost
704 705 706
		size_t l_cnt = 0;
		size_t nl_cnt = 0;
		size_t n_out = 0;
707 708

		auto it2 = vd.getIterator();
Pietro Incardona's avatar
Pietro Incardona committed
709
		count_local_n_local<3,vector_dist<3,float, Point_test<float> >>(vd,it2,bc,box,dom_ext,l_cnt,nl_cnt,n_out);
710 711 712 713 714 715

		// No particles should be out of domain + ghost
		BOOST_REQUIRE_EQUAL(n_out,0ul);

		// Ghost must populated because we synchronized them
		if (k > 524288)
716
		{
717
			BOOST_REQUIRE(nl_cnt != 0);
718 719
			BOOST_REQUIRE(l_cnt > nl_cnt);
		}
720 721 722 723 724 725 726 727 728 729 730

		// Sum all the particles inside the domain
		v_cl.sum(l_cnt);
		v_cl.execute();
		BOOST_REQUIRE_EQUAL(l_cnt,(size_t)k);

		l_cnt = 0;
		nl_cnt = 0;

		// Iterate only on the ghost particles
		auto itg = vd.getGhostIterator();
Pietro Incardona's avatar
Pietro Incardona committed
731
		count_local_n_local<3,vector_dist<3,float, Point_test<float> > >(vd,itg,bc,box,dom_ext,l_cnt,nl_cnt,n_out);
732 733 734 735 736 737

		// No particle on the ghost must be inside the domain
		BOOST_REQUIRE_EQUAL(l_cnt,0ul);

		// Ghost must be populated
		if (k > 524288)
738
		{
739
			BOOST_REQUIRE(nl_cnt != 0);
740
		}
741 742 743 744 745
	}
}

BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_random_walk )
{
746
	Vcluster & v_cl = create_vcluster();
747 748 749 750 751 752

    // 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);
753
	
754 755
	size_t nsz[] = {0,32,4};
	nsz[0] = 65536 * v_cl.getProcessingUnits();
756

Pietro Incardona's avatar
Pietro Incardona committed
757
	print_test_v( "Testing 3D random walk vector k<=",nsz[0]);
758 759

	// 3D test
760
	for (size_t i = 0 ; i < 3 ; i++ )
761
	{
Pietro Incardona's avatar
Pietro Incardona committed
762
		size_t k = nsz[i];
763

764 765 766 767 768 769 770
		BOOST_TEST_CHECKPOINT( "Testing 3D random walk vector 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};

771
		// factor
772
		float factor = pow(create_vcluster().getProcessingUnits()/2.0f,1.0f/3.0f);
773

774
		// ghost
775
		Ghost<3,float> ghost(0.01 / factor);
776 777

		// Distributed vector
incardon's avatar
incardon committed
778
		vector_dist<3,float, Point_test<float> > vd(k,box,bc,ghost);
779 780

		auto it = vd.getIterator();
781 782 783 784 785

		while (it.isNext())
		{
			auto key = it.get();

Pietro Incardona's avatar
Pietro Incardona committed
786 787 788
			vd.getPos(key)[0] = ud(eg);
			vd.getPos(key)[1] = ud(eg);
			vd.getPos(key)[2] = ud(eg);
789 790 791 792

			++it;
		}

793
		vd.map();
794

795
		// 10 step random walk
796

797
		for (size_t j = 0 ; j < 4 ; j++)
798 799 800 801 802 803 804
		{
			auto it = vd.getDomainIterator();

			while (it.isNext())
			{
				auto key = it.get();

Pietro Incardona's avatar
Pietro Incardona committed
805 806 807
				vd.getPos(key)[0] += 0.02 * ud(eg);
				vd.getPos(key)[1] += 0.02 * ud(eg);
				vd.getPos(key)[2] += 0.02 * ud(eg);
808 809 810 811 812 813

				++it;
			}

			vd.map();

814
			vd.ghost_get<0>();
815 816 817 818 819 820

			// Count the local particles and check that the total number is consistent
			size_t cnt = total_n_part_lc(vd,bc);

			BOOST_REQUIRE_EQUAL((size_t)k,cnt);
		}
821 822 823
	}
}

824 825 826 827 828 829 830
BOOST_AUTO_TEST_CASE( vector_dist_periodic_map )
{
	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};

831
	// factor
832
	float factor = pow(create_vcluster().getProcessingUnits()/2.0f,1.0f/3.0f);
833

834
	// ghost
835
	Ghost<3,float> ghost(0.05 / factor);
836 837

	// Distributed vector
incardon's avatar
incardon committed
838
	vector_dist<3,float, Point_test<float> > vd(1,box,bc,ghost);
839 840 841 842 843 844 845 846 847

	// put particles al 1.0, check that they go to 0.0

	auto it = vd.getIterator();

	while (it.isNext())
	{
		auto key = it.get();

Pietro Incardona's avatar
Pietro Incardona committed
848 849 850
		vd.getPos(key)[0] = 1.0;
		vd.getPos(key)[1] = 1.0;
		vd.getPos(key)[2] = 1.0;
851 852 853 854 855 856 857 858 859 860 861 862

		++it;
	}

	vd.map();

	auto it2 = vd.getIterator();

	while (it2.isNext())
	{
		auto key = it2.get();

Pietro Incardona's avatar
Pietro Incardona committed
863
		float f = vd.getPos(key)[0];
864
		BOOST_REQUIRE_EQUAL(f, 0.0);
Pietro Incardona's avatar
Pietro Incardona committed
865
		f = vd.getPos(key)[1];
866
		BOOST_REQUIRE_EQUAL(f, 0.0);
Pietro Incardona's avatar
Pietro Incardona committed
867
		f = vd.getPos(key)[2];
868 869 870 871 872 873
		BOOST_REQUIRE_EQUAL(f, 0.0);

		++it2;
	}
}

874

875 876 877 878 879 880 881
BOOST_AUTO_TEST_CASE( vector_dist_not_periodic_map )
{
	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};

882
	// factor
883
	float factor = pow(create_vcluster().getProcessingUnits()/2.0f,1.0f/3.0f);
884

885
	// ghost
886
	Ghost<3,float> ghost(0.05 / factor);
887 888

	// Distributed vector
incardon's avatar
incardon committed
889
	vector_dist<3,float, Point_test<float> > vd(1,box,bc,ghost);
890 891 892 893 894 895 896 897 898

	// put particles al 1.0, check that they go to 0.0

	auto it = vd.getIterator();

	while (it.isNext())
	{
		auto key = it.get();

Pietro Incardona's avatar
Pietro Incardona committed
899 900 901
		vd.getPos(key)[0] = 1.0;
		vd.getPos(key)[1] = 1.0;
		vd.getPos(key)[2] = 1.0;
902 903 904 905 906 907 908 909 910 911 912 913

		++it;
	}

	vd.map();

	auto it2 = vd.getIterator();

	while (it2.isNext())
	{
		auto key = it2.get();

Pietro Incardona's avatar
Pietro Incardona committed
914
		float f = vd.getPos(key)[0];
915
		BOOST_REQUIRE_EQUAL(f, 1.0);
Pietro Incardona's avatar
Pietro Incardona committed
916
		f = vd.getPos(key)[1];
917
		BOOST_REQUIRE_EQUAL(f, 1.0);
Pietro Incardona's avatar
Pietro Incardona committed
918
		f = vd.getPos(key)[2];
919 920 921 922 923 924
		BOOST_REQUIRE_EQUAL(f, 1.0);

		++it2;
	}
}

925 926
BOOST_AUTO_TEST_CASE( vector_dist_out_of_bound_policy )
{
927
	Vcluster & v_cl = create_vcluster();
928 929 930 931

	if (v_cl.getProcessingUnits() > 8)
		return;

Pietro Incardona's avatar
Pietro Incardona committed
932 933 934
	Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});

	// Boundary conditions
935 936
	size_t bc[3]={NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};

Pietro Incardona's avatar
Pietro Incardona committed
937
	// factor
938
	float factor = pow(create_vcluster().getProcessingUnits()/2.0f,1.0f/3.0f);
Pietro Incardona's avatar
Pietro Incardona committed
939 940 941 942 943

	// ghost
	Ghost<3,float> ghost(0.05 / factor);

	// Distributed vector
incardon's avatar
incardon committed
944
	vector_dist<3,float, Point_test<float> > vd(100,box,bc,ghost);
945 946 947 948 949

	// put particles at out of the boundary, they must be detected and and killed

	auto it = vd.getIterator();

950 951
	size_t cnt = 0;

952 953 954 955
	while (it.isNext())
	{
		auto key = it.get();

956 957
		if (cnt < 1)
		{
Pietro Incardona's avatar
Pietro Incardona committed
958 959 960
			vd.getPos(key)[0] = -0.06;
			vd.getPos(key)[1] = -0.06;
			vd.getPos(key)[2] = -0.06;
961 962 963
		}
		else
		{
Pietro Incardona's avatar
Pietro Incardona committed
964 965 966
			vd.getPos(key)[0] = 0.06;
			vd.getPos(key)[1] = 0.06;
			vd.getPos(key)[2] = 0.06;
967
		}
968

969
		cnt++;
970 971 972
		++it;
	}

Pietro Incardona's avatar
Pietro Incardona committed
973 974
	vd.map();

975 976
	// Particles out of the boundary are killed

977
	size_t cnt_l = vd.size_local();
978

979
	v_cl.sum(cnt_l);
980 981
	v_cl.execute();

982 983 984
	BOOST_REQUIRE_EQUAL(cnt_l,100-v_cl.getProcessingUnits());
}

985
void Test_interacting(Box<3,float> & box)
986
{
987
	Vcluster & v_cl = create_vcluster();
988 989 990 991

	if (v_cl.getProcessingUnits() > 8)
		return;

Pietro Incardona's avatar
Pietro Incardona committed
992
    // set the seed
993 994
	// create the random generator engine
	std::srand(v_cl.getProcessUnitID());
Pietro Incardona's avatar
Pietro Incardona committed
995
    std::default_random_engine eg;
996
    std::uniform_real_distribution<float> ud(-0.5f, 0.5f);
997 998

	size_t nsz[] = {0,32,4};
999 1000 1001 1002

#ifdef TEST_COVERAGE_MODE
	nsz[0] = 5536 * v_cl.getProcessingUnits();
#else
Pietro Incardona's avatar
Pietro Incardona committed
1003
	nsz[0] = 65536 * v_cl.getProcessingUnits();
1004
#endif
Pietro Incardona's avatar
Pietro Incardona committed
1005 1006

	print_test_v("Testing 3D random walk interacting particles vector k=", nsz[0]);
1007 1008 1009 1010 1011 1012

	// 3D test
	for (size_t i = 0 ; i < 3 ; i++ )
	{
		size_t k = nsz[i];

Pietro Incardona's avatar
Pietro Incardona committed
1013
		BOOST_TEST_CHECKPOINT( "Testing 3D random walk interacting particles vector k=" << k );
1014 1015 1016 1017 1018

		// Boundary conditions
		size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};

		// factor
1019
		float factor = pow(create_vcluster().getProcessingUnits()/2.0f,1.0f/3.0f);
1020 1021 1022 1023 1024 1025 1026 1027

		// interaction radius
		float r_cut = 0.01 / factor;

		// ghost
		Ghost<3,float> ghost(r_cut);

		// Distributed vector
incardon's avatar
incardon committed
1028
		vector_dist<3,float, Point_test<float> > vd(k,box,bc,ghost);
1029 1030 1031 1032 1033 1034 1035

		auto it = vd.getIterator();

		while (it.isNext())
		{
			auto key = it.get();

Pietro Incardona's avatar
Pietro Incardona committed
1036 1037 1038
			vd.getPos(key)[0] = ud(eg);
			vd.getPos(key)[1] = ud(eg);
			vd.getPos(key)[2] = ud(eg);
1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056

			++it;
		}

		vd.map();

		// 4 step random walk

		for (size_t j = 0 ; j < 4 ; j++)
		{
			auto it = vd.getDomainIterator();

			// Move the particles

			while (it.isNext())
			{
				auto key = it.get();

Pietro Incardona's avatar
Pietro Incardona committed
1057 1058 1059
				vd.getPos(key)[0] += 0.02 * ud(eg);
				vd.getPos(key)[1] += 0.02 * ud(eg);
				vd.getPos(key)[2] += 0.02 * ud(eg);
1060 1061 1062 1063 1064 1065

				++it;
			}

			vd.map();

1066 1067
			vd.write("Without_ghost");

1068 1069
			vd.ghost_get<0>();

1070 1071 1072
			vd.write("With_ghost");
			vd.getDecomposition().write("With_dec_ghost");

1073 1074
			// get the cell list with a cutoff radius

1075 1076
			bool error = false;

1077 1078 1079 1080 1081 1082 1083 1084 1085 1086
			auto NN = vd.getCellList(0.01 / factor);

			// iterate across the domain particle

			auto it2 = vd.getDomainIterator();

			while (it2.isNext())
			{
				auto p = it2.get();

Pietro Incardona's avatar
Pietro Incardona committed
1087
				Point<3,float> xp = vd.getPos(p);
1088

incardon's avatar
incardon committed
1089
				auto Np = NN.getCellIterator(NN.getCell(xp));
1090 1091 1092 1093 1094 1095 1096

				while (Np.isNext())
				{
					auto q = Np.get();

					// repulsive

Pietro Incardona's avatar
Pietro Incardona committed
1097
					Point<3,float> xq = vd.getPos(q);
1098 1099 1100 1101
					Point<3,float> f = (xp - xq);

					float distance = f.norm();

1102
					// Particle should be inside 2 * r_cut range
1103

1104 1105
					if (distance > 2*r_cut*sqrt(2))
						error = true;
1106 1107 1108 1109 1110 1111 1112

					++Np;
				}

				++it2;
			}

1113 1114 1115 1116
			// Error

			BOOST_REQUIRE_EQUAL(error,false);

1117 1118 1119 1120 1121 1122
			// Count the local particles and check that the total number is consistent
			size_t cnt = total_n_part_lc(vd,bc);

			BOOST_REQUIRE_EQUAL((size_t)k,cnt);
		}
	}
1123
}
1124

1125 1126 1127 1128 1129 1130 1131 1132 1133
BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_interacting_particles )
{
	Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
	Test_interacting(box);

	Box<3,float> box2({-0.5,-0.5,-0.5},{0.5,0.5,0.5});
	Test_interacting(box2);
}

1134 1135
BOOST_AUTO_TEST_CASE( vector_dist_grid_iterator )
{
1136 1137 1138
#ifdef TEST_COVERAGE_MODE
	long int k = 32*32*32*create_vcluster().getProcessingUnits();
#else
1139
	long int k = 64*64*64*create_vcluster().getProcessingUnits();
1140
#endif
1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167
	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
incardon's avatar
incardon committed
1168
		vector_dist<3,float, Point_test<float> > vd(0,box,bc,ghost);
1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202

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

1203 1204
BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
{
1205 1206 1207
#ifdef TEST_COVERAGE_MODE
	long int k = 16*16*16*create_vcluster().getProcessingUnits();
#else
1208
	long int k = 64*64*64*create_vcluster().getProcessingUnits();
1209
#endif
1210
	k = std::pow(k, 1/3.);
1211

1212 1213 1214
	long int big_step = k / 30;
	big_step = (big_step == 0)?1:big_step;
	long int small_step = 21;
1215

1216
	print_test( "Testing cell and verlet list k<=",k);
1217

1218
	// 3D test
1219
	for ( ; k > 8*big_step ; k-= (k > 2*big_step)?big_step:small_step )
1220
	{
1221
		Vcluster & v_cl = create_vcluster();
1222

1223
		const size_t Ng = k;
1224

1225 1226
		// we create a 128x128x128 Grid iterator
		size_t sz[3] = {Ng,Ng,Ng};
1227

1228
		Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
1229

1230 1231
		// Boundary conditions
		size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
1232

Pietro Incardona's avatar
Pietro Incardona committed
1233 1234 1235 1236 1237
		float spacing = 1.0/Ng;
		float first_dist = spacing;
		float second_dist = sqrt(2.0*spacing*spacing);
		float third_dist = sqrt(3.0 * spacing*spacing);

1238
		// ghost
Pietro Incardona's avatar
Pietro Incardona committed
1239
		Ghost<3,float> ghost(third_dist*1.1);
1240

1241
		// Distributed vector
incardon's avatar
incardon committed
1242
		vector_dist<3,float, Point_test<float> > vd(0,box,bc,ghost);
1243

1244 1245
		// Put particles on a grid creating a Grid iterator
		auto it = vd.getGridIterator(sz);
1246

1247 1248 1249
		while (it.isNext())
		{
			vd.add();
1250

1251
			auto key = it.get();
1252

Pietro Incardona's avatar
Pietro Incardona committed
1253 1254 1255
			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);
1256

1257 1258
			++it;
		}
1259

1260 1261 1262 1263
		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);

1264 1265
		// distribute particles and sync ghost
		vd.map();
1266

1267 1268 1269 1270
		// 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();
1271

1272
		BOOST_REQUIRE_EQUAL(total,(Ng) * (Ng) * (Ng));
1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286

		vd.ghost_get<0>();

		// calculate the distance of the first, second and third neighborhood particle
		// Consider that they are on a regular grid

		// add a 5% to dist

		first_dist += first_dist * 0.05;
		second_dist += second_dist * 0.05;
		third_dist += third_dist * 0.05;

		// Create a verlet list for each particle

1287
		VerletList<3,float,FAST,shift<3,float>> verlet = vd.getVerlet(third_dist);
1288

1289
		bool correct = true;
1290

Pietro Incardona's avatar
Pietro Incardona committed
1291 1292
		BOOST_REQUIRE_EQUAL(vd.size_local(),verlet.size());

1293 1294
		// for each particle
		for (size_t i = 0 ; i < verlet.size() ; i++)
1295
		{
1296 1297 1298 1299
			// first NN
			size_t first_NN = 0;
			size_t second_NN = 0;
			size_t third_NN = 0;
1300

Pietro Incardona's avatar
Pietro Incardona committed
1301
			Point<3,float> p = vd.getPos(i);
1302

1303
			// for each neighborhood particle
1304
			for (size_t j = 0 ; j < verlet.getNNPart(i) ; j++)
1305
			{
1306
				Point<3,float> q = vd.getPos(verlet.get(i,j));
1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317

				float dist = p.distance(q);

				if (dist <= first_dist)
					first_NN++;
				else if (dist <= second_dist)
					second_NN++;
				else
					third_NN++;
			}