vector_dist_unit_test.hpp 25.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_

incardon's avatar
incardon committed
11
#include <random>
incardon's avatar
incardon committed
12 13
#include "Vector/vector_dist.hpp"

14 15 16 17 18 19 20 21 22 23 24 25 26 27
/*! \brief Count the total number of particles
 *
 * \param vd distributed vector
 * \param bc boundary conditions
 *
 */
template<unsigned int dim> size_t total_n_part_lc(vector_dist<dim,float, Point_test<float>, CartDecomposition<dim,float> > & vd, size_t (& bc)[dim])
{
	typedef Point<dim,float> s;

	Vcluster & v_cl = vd.getVC();
	auto it2 = vd.getDomainIterator();
	const CartDecomposition<3,float> & ct = vd.getDecomposition();

28 29
	bool noOut = true;

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

35
		noOut &= ct.isLocal(vd.template getPos<s::x>(key));
36 37 38 39 40 41

		cnt++;

		++it2;
	}

42 43
	BOOST_REQUIRE_EQUAL(noOut,true);

44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
	//
	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
 *
 */
template<unsigned int dim> inline void count_local_n_local(vector_dist<dim,float, Point_test<float>, CartDecomposition<dim,float> > & 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)
{
	typedef Point<dim,float> s;
	const CartDecomposition<dim,float> & ct = vd.getDecomposition();

	while (it.isNext())
	{
		auto key = it.get();
		// Check if it is in the domain
		if (box.isInsideNP(vd.template getPos<s::x>(key)) == true)
		{
			// Check if local
			if (ct.isLocalBC(vd.template getPos<s::x>(key),bc) == true)
				l_cnt++;
			else
				nl_cnt++;
		}
		else
		{
			nl_cnt++;
		}

		// Check that all particles are inside the Domain + Ghost part
		if (dom_ext.isInside(vd.template getPos<s::x>(key)) == false)
				n_out++;

		++it;
	}
}

incardon's avatar
incardon committed
93 94
BOOST_AUTO_TEST_SUITE( vector_dist_test )

95 96 97 98 99 100
void print_test(std::string test, size_t sz)
{
	if (global_v_cluster->getProcessUnitID() == 0)
		std::cout << test << " " << sz << "\n";
}

101 102 103 104 105 106 107 108
BOOST_AUTO_TEST_CASE( vector_dist_ghost )
{
	// Communication object
	Vcluster & v_cl = *global_v_cluster;

	typedef Point_test<float> p;
	typedef Point<2,float> s;

109 110 111 112
	// Get the default minimum number of sub-sub-domain per processor (granularity of the decomposition)
	size_t n_sub = vector_dist<2,float, Point_test<float>, CartDecomposition<2,float> >::getDefaultNsubsub() * v_cl.getProcessingUnits();
	// 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
113

114 115
	//! [Create a vector of elements distributed on a grid like way]

116
	Box<2,float> box({0.0,0.0},{1.0,1.0});
117
	size_t g_div[]= {sz,sz};
118

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

122
	// Calculate the number of elements this processor is going to obtain
123 124 125 126 127 128 129 130
	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);

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

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

	// Calculate the grid spacing
	Point<2,float> spacing = box.getP2();
	spacing = spacing / g_div;

	// middle spacing
	Point<2,float> m_spacing = spacing / 2;

incardon's avatar
incardon committed
145 146 147
	// 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);

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

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

	// size_t
	size_t cobj = 0;
156

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

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

		// set the particle position

		vd.template getPos<s::x>(key_v)[0] = key.get(0) * spacing[0] + m_spacing[0];
		vd.template getPos<s::x>(key_v)[1] = key.get(1) * spacing[1] + m_spacing[1];

170 171 172
		cobj++;

		++v_it;
173 174 175
		++it;
	}

176 177 178
	//! [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
179 180 181 182 183
	// predicted one
	BOOST_REQUIRE_EQUAL(v_it.isNext(),false);
	BOOST_REQUIRE_EQUAL(it.isNext(),false);
	BOOST_REQUIRE_EQUAL(cobj,p_np);

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

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

189
	auto v_it2 = vd.getIterator();
190

191
	while (v_it2.isNext())
192
	{
193
		auto key = v_it2.get();
194 195 196 197 198 199 200

		// fill with the processor ID where these particle live
		vd.template getProp<p::s>(key) = vd.getPos<s::x>(key)[0] + vd.getPos<s::x>(key)[1] * 16;
		vd.template getProp<p::v>(key)[0] = v_cl.getProcessUnitID();
		vd.template getProp<p::v>(key)[1] = v_cl.getProcessUnitID();
		vd.template getProp<p::v>(key)[2] = v_cl.getProcessUnitID();

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

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

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

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

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

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

218 219
	size_t n_part = 0;

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

		// Check the received data
		BOOST_REQUIRE_EQUAL(vd.getPos<s::x>(key)[0] + vd.getPos<s::x>(key)[1] * 16,vd.template getProp<p::s>(key));

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

232
		// check if the received data are in one of the ghost boxes
233
		for ( ; b < dec.getNEGhostBox() ; b++)
234
		{
235
			if (dec.getEGhostBox(b).isInside(vd.getPos<s::x>(key)) == true )
236 237 238 239 240 241 242
			{
				is_in = true;

				// Add
				vb.get(b)++;
				lb = b;
			}
243 244 245 246
		}
		BOOST_REQUIRE_EQUAL(is_in,true);

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

249
		n_part++;
250 251 252
		++g_it;
	}

253 254 255 256
	if (v_cl.getProcessingUnits() > 1)
	{
		BOOST_REQUIRE(n_part != 0);
	}
257

258 259 260 261 262
    CellDecomposer_sm<2,float> cd(SpaceBox<2,float>(box),g_div,0);

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

265 266 267 268 269
		if (n_point != vb.get(i))
		{
			std::cout << n_point << "  " << dec.getEGhostBoxProcessor(i) << "  " << v_cl.getProcessUnitID() << dec.getEGhostBox(i).toString() << "\n";
		}
		//BOOST_REQUIRE_EQUAL(n_point,vb.get(i));
270
	}
incardon's avatar
incardon committed
271 272
}

273 274 275 276 277 278
void print_test_v(std::string test, size_t sz)
{
	if (global_v_cluster->getProcessUnitID() == 0)
		std::cout << test << " " << sz << "\n";
}

incardon's avatar
incardon committed
279 280 281 282 283 284 285 286 287 288 289 290 291 292
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;
}

293
BOOST_AUTO_TEST_CASE( vector_dist_iterator_test_use_2d )
incardon's avatar
incardon committed
294 295 296 297 298 299 300 301 302 303 304
{
	typedef Point<2,float> s;

	Vcluster & v_cl = *global_v_cluster;

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

305
    long int k = 524288 * v_cl.getProcessingUnits();
incardon's avatar
incardon committed
306

307
	long int big_step = k / 4;
308
	big_step = (big_step == 0)?1:big_step;
incardon's avatar
incardon committed
309

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

312
	// 2D test
incardon's avatar
incardon committed
313
	for ( ; k >= 2 ; k-= decrement(k,big_step) )
incardon's avatar
incardon committed
314
	{
315
		BOOST_TEST_CHECKPOINT( "Testing 2D vector k=" << k );
316 317 318

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

319
		Box<2,float> box({0.0,0.0},{1.0,1.0});
320 321 322 323 324

		// 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));
incardon's avatar
incardon committed
325

326
		auto it = vd.getIterator();
incardon's avatar
incardon committed
327

328 329 330 331 332 333 334 335 336 337 338 339
		while (it.isNext())
		{
			auto key = it.get();

			vd.template getPos<s::x>(key)[0] = ud(eg);
			vd.template getPos<s::x>(key)[1] = ud(eg);

			++it;
		}

		vd.map();

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

342 343 344
		// Check if we have all the local particles
		size_t cnt = 0;
		const CartDecomposition<2,float> & ct = vd.getDecomposition();
345
		auto it2 = vd.getIterator();
346

347
		while (it2.isNext())
348
		{
349
			auto key = it2.get();
350 351 352 353 354 355

			// Check if local
			BOOST_REQUIRE_EQUAL(ct.isLocal(vd.template getPos<s::x>(key)),true);

			cnt++;

356
			++it2;
357 358 359 360 361
		}

		//
		v_cl.sum(cnt);
		v_cl.execute();
362
		BOOST_REQUIRE_EQUAL((long int)cnt,k);
incardon's avatar
incardon committed
363
	}
364
}
incardon's avatar
incardon committed
365

366 367 368
BOOST_AUTO_TEST_CASE( vector_dist_iterator_test_use_3d )
{
	typedef Point<3,float> s;
incardon's avatar
incardon committed
369

370
	Vcluster & v_cl = *global_v_cluster;
incardon's avatar
incardon committed
371

372 373 374 375 376 377
    // 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);

378
    long int k = 524288 * v_cl.getProcessingUnits();
379

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

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

385
	// 3D test
incardon's avatar
incardon committed
386
	for ( ; k >= 2 ; k-= decrement(k,big_step) )
incardon's avatar
incardon committed
387
	{
388
		BOOST_TEST_CHECKPOINT( "Testing 3D vector k=" << k );
389 390 391

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

392
		Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
393 394 395 396 397

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

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

399
		auto it = vd.getIterator();
incardon's avatar
incardon committed
400

401 402 403
		while (it.isNext())
		{
			auto key = it.get();
incardon's avatar
incardon committed
404

405 406 407 408 409 410 411 412 413
			vd.template getPos<s::x>(key)[0] = ud(eg);
			vd.template getPos<s::x>(key)[1] = ud(eg);
			vd.template getPos<s::x>(key)[2] = ud(eg);

			++it;
		}

		vd.map();

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

416 417 418
		// Check if we have all the local particles
		size_t cnt = 0;
		const CartDecomposition<3,float> & ct = vd.getDecomposition();
419
		auto it2 = vd.getIterator();
420

421
		while (it2.isNext())
422
		{
423
			auto key = it2.get();
424 425

			// Check if local
incardon's avatar
incardon committed
426
			BOOST_REQUIRE_EQUAL(ct.isLocal(vd.template getPos<s::x>(key)),true);
427 428 429

			cnt++;

430
			++it2;
431 432 433 434 435
		}

		//
		v_cl.sum(cnt);
		v_cl.execute();
436
		BOOST_REQUIRE_EQUAL(cnt,(size_t)k);
437
	}
incardon's avatar
incardon committed
438 439
}

440 441 442 443 444 445 446 447 448 449 450 451 452 453
BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_use_2d )
{
	typedef Point<2,float> s;

	Vcluster & v_cl = *global_v_cluster;

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

454
	long int big_step = k / 4;
455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470
	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 );

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

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

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

471 472 473
		// factor
		float factor = pow(global_v_cluster->getProcessingUnits()/2.0f,1.0f/3.0f);

474
		// ghost
475
		Ghost<2,float> ghost(0.01 / factor);
476

477
		// ghost2 (a little bigger because of round off error)
478
		Ghost<2,float> ghost2(0.05001 / factor);
479 480

		// Distributed vector
481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496
		vector_dist<2,float, Point_test<float>, CartDecomposition<2,float> > vd(k,box,bc,ghost);

		auto it = vd.getIterator();

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

			vd.template getPos<s::x>(key)[0] = ud(eg);
			vd.template getPos<s::x>(key)[1] = ud(eg);

			++it;
		}

		vd.map();

497
		// sync the ghost, only the property zero
498 499 500 501
		vd.ghost_get<0>();

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

502
		// Domain + ghost box
503
		Box<2,float> dom_ext = box;
504
		dom_ext.enlarge(ghost2);
505

506 507 508 509
		// Iterate on all particles domain + ghost
		size_t l_cnt = 0;
		size_t nl_cnt = 0;
		size_t n_out = 0;
510 511


512 513
		auto it2 = vd.getIterator();
		count_local_n_local(vd,it2,bc,box,dom_ext,l_cnt,nl_cnt,n_out);
514

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

518
		// Ghost must populated because we synchronized them
519
		if (k > 524288)
520
		{
521
			BOOST_REQUIRE(nl_cnt != 0);
522 523
			BOOST_REQUIRE(l_cnt > nl_cnt);
		}
524

525
		// Sum all the particles inside the domain
526 527
		v_cl.sum(l_cnt);
		v_cl.execute();
528 529

		// count that they are equal to the initial total number
530
		BOOST_REQUIRE_EQUAL((long int)l_cnt,k);
531 532 533 534 535 536 537 538 539 540 541 542 543

		l_cnt = 0;
		nl_cnt = 0;

		// Iterate only on the ghost particles
		auto itg = vd.getGhostIterator();
		count_local_n_local(vd,itg,bc,box,dom_ext,l_cnt,nl_cnt,n_out);

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

		// Ghost must be populated
		if (k > 524288)
544
		{
545
			BOOST_REQUIRE(nl_cnt != 0);
546
		}
547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563
	}
}

BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_use_3d )
{
	typedef Point<3,float> s;

	Vcluster & v_cl = *global_v_cluster;

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

564
	long int big_step = k / 4;
565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580
	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 );

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

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

581 582 583
		// factor
		float factor = pow(global_v_cluster->getProcessingUnits()/2.0f,1.0f/3.0f);

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

587
		// ghost2 (a little bigger because of round off error)
588
		Ghost<3,float> ghost2(0.05001 / factor);
589 590

		// Distributed vector
591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612
		vector_dist<3,float, Point_test<float>, CartDecomposition<3,float> > vd(k,box,bc,ghost);

		auto it = vd.getIterator();

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

			vd.template getPos<s::x>(key)[0] = ud(eg);
			vd.template getPos<s::x>(key)[1] = ud(eg);
			vd.template getPos<s::x>(key)[2] = ud(eg);

			++it;
		}

		vd.map();

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

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

613 614 615 616 617
		// Domain + ghost
		Box<3,float> dom_ext = box;
		dom_ext.enlarge(ghost2);

		// Iterate on all particles domain + ghost
618 619 620
		size_t l_cnt = 0;
		size_t nl_cnt = 0;
		size_t n_out = 0;
621 622 623 624 625 626 627 628 629

		auto it2 = vd.getIterator();
		count_local_n_local(vd,it2,bc,box,dom_ext,l_cnt,nl_cnt,n_out);

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

		// Ghost must populated because we synchronized them
		if (k > 524288)
630
		{
631
			BOOST_REQUIRE(nl_cnt != 0);
632 633
			BOOST_REQUIRE(l_cnt > nl_cnt);
		}
634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651

		// 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();
		count_local_n_local(vd,itg,bc,box,dom_ext,l_cnt,nl_cnt,n_out);

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

		// Ghost must be populated
		if (k > 524288)
652
		{
653
			BOOST_REQUIRE(nl_cnt != 0);
654
		}
655 656 657 658 659 660 661 662 663 664 665 666 667 668
	}
}

BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_random_walk )
{
	typedef Point<3,float> s;

	Vcluster & v_cl = *global_v_cluster;

    // 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);
669
	
670 671
	size_t nsz[] = {0,32,4};
	nsz[0] = 65536 * v_cl.getProcessingUnits();
672

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

	// 3D test
676
	for (size_t i = 0 ; i < 3 ; i++ )
677
	{
Pietro Incardona's avatar
Pietro Incardona committed
678
		size_t k = nsz[i];
679

680 681 682 683 684 685 686
		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};

687 688 689
		// factor
		float factor = pow(global_v_cluster->getProcessingUnits()/2.0f,1.0f/3.0f);

690
		// ghost
691
		Ghost<3,float> ghost(0.01 / factor);
692 693 694 695 696

		// Distributed vector
		vector_dist<3,float, Point_test<float>, CartDecomposition<3,float> > vd(k,box,bc,ghost);

		auto it = vd.getIterator();
697 698 699 700 701

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

702 703 704
			vd.template getPos<s::x>(key)[0] = ud(eg);
			vd.template getPos<s::x>(key)[1] = ud(eg);
			vd.template getPos<s::x>(key)[2] = ud(eg);
705 706 707 708

			++it;
		}

709
		vd.map();
710

711
		// 10 step random walk
712

713
		for (size_t j = 0 ; j < 4 ; j++)
714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729
		{
			auto it = vd.getDomainIterator();

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

				vd.template getPos<s::x>(key)[0] += 0.02 * ud(eg);
				vd.template getPos<s::x>(key)[1] += 0.02 * ud(eg);
				vd.template getPos<s::x>(key)[2] += 0.02 * ud(eg);

				++it;
			}

			vd.map();

730
			vd.ghost_get<0>();
731 732 733 734 735 736

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

740 741 742 743 744 745 746 747 748
BOOST_AUTO_TEST_CASE( vector_dist_periodic_map )
{
	typedef Point<3,float> s;

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

749 750 751
	// factor
	float factor = pow(global_v_cluster->getProcessingUnits()/2.0f,1.0f/3.0f);

752
	// ghost
753
	Ghost<3,float> ghost(0.05 / factor);
754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800

	// Distributed vector
	vector_dist<3,float, Point_test<float>, CartDecomposition<3,float> > vd(1,box,bc,ghost);

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

	auto it = vd.getIterator();

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

		vd.template getPos<s::x>(key)[0] = 1.0;
		vd.template getPos<s::x>(key)[1] = 1.0;
		vd.template getPos<s::x>(key)[2] = 1.0;

		++it;
	}

	vd.map();

	auto it2 = vd.getIterator();

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

		float f = vd.template getPos<s::x>(key)[0];
		BOOST_REQUIRE_EQUAL(f, 0.0);
		f = vd.template getPos<s::x>(key)[1];
		BOOST_REQUIRE_EQUAL(f, 0.0);
		f = vd.template getPos<s::x>(key)[2];
		BOOST_REQUIRE_EQUAL(f, 0.0);

		++it2;
	}
}

BOOST_AUTO_TEST_CASE( vector_dist_not_periodic_map )
{
	typedef Point<3,float> s;

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

801 802 803
	// factor
	float factor = pow(global_v_cluster->getProcessingUnits()/2.0f,1.0f/3.0f);

804
	// ghost
805
	Ghost<3,float> ghost(0.05 / factor);
806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843

	// Distributed vector
	vector_dist<3,float, Point_test<float>, CartDecomposition<3,float> > vd(1,box,bc,ghost);

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

	auto it = vd.getIterator();

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

		vd.template getPos<s::x>(key)[0] = 1.0;
		vd.template getPos<s::x>(key)[1] = 1.0;
		vd.template getPos<s::x>(key)[2] = 1.0;

		++it;
	}

	vd.map();

	auto it2 = vd.getIterator();

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

		float f = vd.template getPos<s::x>(key)[0];
		BOOST_REQUIRE_EQUAL(f, 1.0);
		f = vd.template getPos<s::x>(key)[1];
		BOOST_REQUIRE_EQUAL(f, 1.0);
		f = vd.template getPos<s::x>(key)[2];
		BOOST_REQUIRE_EQUAL(f, 1.0);

		++it2;
	}
}

844 845
BOOST_AUTO_TEST_CASE( vector_dist_out_of_bound_policy )
{
846
	Vcluster & v_cl = *global_v_cluster;
847 848 849 850 851 852

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

	typedef Point<3,float> s;

Pietro Incardona's avatar
Pietro Incardona committed
853 854 855
	Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});

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

Pietro Incardona's avatar
Pietro Incardona committed
858 859 860 861 862 863 864
	// factor
	float factor = pow(global_v_cluster->getProcessingUnits()/2.0f,1.0f/3.0f);

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

	// Distributed vector
865 866 867 868 869 870
	vector_dist<3,float, Point_test<float>, CartDecomposition<3,float> > vd(100,box,bc,ghost);

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

	auto it = vd.getIterator();

871 872
	size_t cnt = 0;

873 874 875 876
	while (it.isNext())
	{
		auto key = it.get();

877 878 879 880 881 882 883 884 885 886 887 888
		if (cnt < 1)
		{
			vd.template getPos<s::x>(key)[0] = -0.06;
			vd.template getPos<s::x>(key)[1] = -0.06;
			vd.template getPos<s::x>(key)[2] = -0.06;
		}
		else
		{
			vd.template getPos<s::x>(key)[0] = 0.06;
			vd.template getPos<s::x>(key)[1] = 0.06;
			vd.template getPos<s::x>(key)[2] = 0.06;
		}
889

890
		cnt++;
891 892 893
		++it;
	}

Pietro Incardona's avatar
Pietro Incardona committed
894 895
	vd.map();

896 897
	// Particles out of the boundary are killed

898
	size_t cnt_l = vd.size_local();
899

900
	v_cl.sum(cnt_l);
901 902
	v_cl.execute();

903 904 905 906 907
	BOOST_REQUIRE_EQUAL(cnt_l,100-v_cl.getProcessingUnits());
}

BOOST_AUTO_TEST_CASE( vector_dist_periodic_test_interacting_particles )
{
Pietro Incardona's avatar
Pietro Incardona committed
908

909 910 911 912 913 914 915
	typedef Point<3,float> s;

	Vcluster & v_cl = *global_v_cluster;

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

Pietro Incardona's avatar
Pietro Incardona committed
916
    // set the seed
917 918
	// create the random generator engine
	std::srand(v_cl.getProcessUnitID());
Pietro Incardona's avatar
Pietro Incardona committed
919 920
    std::default_random_engine eg;
    std::uniform_real_distribution<float> ud(0.0f, 1.0f);
921 922

	size_t nsz[] = {0,32,4};
Pietro Incardona's avatar
Pietro Incardona committed
923 924 925
	nsz[0] = 65536 * v_cl.getProcessingUnits();

	print_test_v("Testing 3D random walk interacting particles vector k=", nsz[0]);
926 927 928 929 930 931

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

Pietro Incardona's avatar
Pietro Incardona committed
932
		BOOST_TEST_CHECKPOINT( "Testing 3D random walk interacting particles vector k=" << k );
933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990

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

		// factor
		float factor = pow(global_v_cluster->getProcessingUnits()/2.0f,1.0f/3.0f);

		// interaction radius
		float r_cut = 0.01 / factor;

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

		// Distributed vector
		vector_dist<3,float, Point_test<float>, CartDecomposition<3,float> > vd(k,box,bc,ghost);

		auto it = vd.getIterator();

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

			vd.template getPos<s::x>(key)[0] = ud(eg);
			vd.template getPos<s::x>(key)[1] = ud(eg);
			vd.template getPos<s::x>(key)[2] = ud(eg);

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

				vd.template getPos<s::x>(key)[0] += 0.02 * ud(eg);
				vd.template getPos<s::x>(key)[1] += 0.02 * ud(eg);
				vd.template getPos<s::x>(key)[2] += 0.02 * ud(eg);

				++it;
			}

			vd.map();

			vd.ghost_get<0>();

			// get the cell list with a cutoff radius

991 992
			bool error = false;

993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017
			auto NN = vd.getCellList(0.01 / factor);

			// iterate across the domain particle

			auto it2 = vd.getDomainIterator();

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

				Point<3,float> xp = vd.getPos<0>(p);

				auto Np = NN.getIterator(NN.getCell(vd.getPos<0>(p)));

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

					// repulsive

					Point<3,float> xq = vd.getPos<0>(q);
					Point<3,float> f = (xp - xq);

					float distance = f.norm();

1018
					// Particle should be inside 2 * r_cut range
1019

1020 1021
					if (distance > 2*r_cut*sqrt(2))
						error = true;
1022 1023 1024 1025 1026 1027 1028

					++Np;
				}

				++it2;
			}

1029 1030 1031 1032
			// Error

			BOOST_REQUIRE_EQUAL(error,false);

1033 1034 1035 1036 1037 1038
			// 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);
		}
	}
1039
}
1040 1041 1042

BOOST_AUTO_TEST_CASE( vector_dist_cell_verlet_test )
{
1043
	long int k = 64*64*64*global_v_cluster->getProcessingUnits();
1044
	k = std::pow(k, 1/3.);
1045

1046 1047 1048
	long int big_step = k / 30;
	big_step = (big_step == 0)?1:big_step;
	long int small_step = 21;
1049

1050
	print_test( "Testing cell and verlet list k<=",k);
1051

1052
	// 3D test
1053
	for ( ; k > 8*big_step ; k-= (k > 2*big_step)?big_step:small_step )
1054 1055
	{
		typedef Point<3,float> s;
1056

1057
		Vcluster & v_cl = *global_v_cluster;
1058

1059
		const size_t Ng = k;
1060

1061 1062
		// we create a 128x128x128 Grid iterator
		size_t sz[3] = {Ng,Ng,Ng};
1063

1064
		Box<3,float> box({0.0,0.0,0.0},{1.0,1.0,1.0});
1065

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

1069 1070
		// ghost
		Ghost<3,float> ghost(1.0/(Ng-2));
1071

1072 1073
		// Distributed vector
		vector_dist<3,float, Point_test<float>, CartDecomposition<3,float> > vd(0,box,bc,ghost);
1074

1075 1076
		// Put particles on a grid creating a Grid iterator
		auto it = vd.getGridIterator(sz);
1077

1078 1079 1080
		while (it.isNext())
		{
			vd.add();
1081

1082
			auto key = it.get();
1083

1084 1085 1086
			vd.template getLastPos<s::x>()[0] = key.get(0) * it.getSpacing(0);
			vd.template getLastPos<s::x>()[1] = key.get(1) * it.getSpacing(1);
			vd.template getLastPos<s::x>()[2] = key.get(2) * it.getSpacing(2);
1087

1088 1089
			++it;
		}
1090

1091 1092
		// distribute particles and sync ghost
		vd.map();
1093

1094 1095 1096 1097
		// 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();
1098

1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119
		BOOST_REQUIRE_EQUAL(total,(Ng-1) * (Ng-1) * (Ng-1));

		vd.ghost_get<0>();

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

		float spacing = it.getSpacing(0);
		float first_dist = spacing;
		float second_dist = sqrt(2.0*spacing*spacing);
		float third_dist = sqrt(3.0 * spacing*spacing);

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

		openfpm::vector<openfpm::vector<size_t>> verlet;
1120

1121
		vd.getVerlet(verlet,third_dist);
1122

1123
		bool correct = true;
1124

1125 1126
		// for each particle
		for (size_t i = 0 ; i < verlet.size() ; i++)
1127
		{
1128 1129 1130 1131
			// first NN
			size_t first_NN = 0;
			size_t second_NN = 0;
			size_t third_NN = 0;
1132

1133
			Point<3,float> p = vd.getPos<0>(i);
1134

1135 1136 1137 1138
			// for each neighborhood particle
			for (size_t j = 0 ; j < verlet.get(i).size() ; j++)
			{
				auto & NN = verlet.get(i);
1139

1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154
				Point<3,float> q = vd.getPos<0>(NN.get(j));

				float dist = p.distance(q);

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

			correct &= (first_NN == 6);
			correct &= (second_NN == 12);
			correct &= (third_NN == 8);
1155 1156
		}

1157
		BOOST_REQUIRE_EQUAL(correct,true);
1158 1159 1160
	}
}

incardon's avatar
incardon committed
1161 1162 1163
BOOST_AUTO_TEST_SUITE_END()

#endif /* VECTOR_DIST_UNIT_TEST_HPP_ */