1/*
2 * vector_dist_cell_list_tests.hpp
3 *
4 * Created on: Aug 16, 2016
5 * Author: i-bird
6 */
7
8#include "config.h"
9#define BOOST_TEST_DYN_LINK
10#include <boost/test/unit_test.hpp>
11#include "Point_test.hpp"
12#include "Vector/performance/vector_dist_performance_common.hpp"
13#include "Vector/vector_dist.hpp"
14
15extern void print_test_v(std::string test, size_t sz);
16extern long int decrement(long int k, long int step);
17
18///////////////////////// test hilb ///////////////////////////////
19
20void test_reorder_sfc(reorder_opt opt)
21{
22 Vcluster<> & v_cl = create_vcluster();
23
24 if (v_cl.getProcessingUnits() > 48)
25 return;
26
27 // set the seed
28 // create the random generator engine
29 std::srand(v_cl.getProcessUnitID());
30 std::default_random_engine eg;
31 std::uniform_real_distribution<float> ud(0.0f, 1.0f);
32
33#ifdef TEST_COVERAGE_MODE
34 long int k = 24288 * v_cl.getProcessingUnits();
35#else
36 long int k = 524288 * v_cl.getProcessingUnits();
37#endif
38
39 long int big_step = k / 4;
40 big_step = (big_step == 0)?1:big_step;
41
42 print_test_v( "Testing 2D vector with sfc curve reordering k<=",k);
43
44 // 2D test
45 for ( ; k >= 2 ; k-= decrement(k,big_step) )
46 {
47 BOOST_TEST_CHECKPOINT( "Testing 2D vector with sfc curve reordering k=" << k );
48
49 Box<2,float> box({0.0,0.0},{1.0,1.0});
50
51 // Boundary conditions
52 size_t bc[2]={NON_PERIODIC,NON_PERIODIC};
53
54 vector_dist<2,float, Point_test<float> > vd(k,box,bc,Ghost<2,float>(0.01));
55
56 auto it = vd.getIterator();
57
58 while (it.isNext())
59 {
60 auto key = it.get();
61
62 vd.getPos(key)[0] = ud(eg);
63 vd.getPos(key)[1] = ud(eg);
64
65 ++it;
66 }
67
68 vd.map();
69
70 // in case of SE_CLASS3 get a cell-list without ghost get is an error
71
72 // Create first cell list
73
74 auto NN1 = vd.getCellList(0.01,true);
75
76 //An order of a curve
77 int32_t m = 6;
78
79 //Reorder a vector
80 vd.reorder(m,opt);
81
82 // Create second cell list
83 auto NN2 = vd.getCellList(0.01,true);
84
85 //Check equality of cell sizes
86 for (size_t i = 0 ; i < NN1.getGrid().size() ; i++)
87 {
88 size_t n1 = NN1.getNelements(i);
89 size_t n2 = NN2.getNelements(i);
90
91 BOOST_REQUIRE_EQUAL(n1,n2);
92 }
93 }
94}
95
96///////////////////////// test hilb ///////////////////////////////
97
98void test_reorder_cl()
99{
100 Vcluster<> & v_cl = create_vcluster();
101
102 if (v_cl.getProcessingUnits() > 48)
103 return;
104
105 // set the seed
106 // create the random generator engine
107 std::srand(v_cl.getProcessUnitID());
108 std::default_random_engine eg;
109 std::uniform_real_distribution<float> ud(0.0f, 1.0f);
110
111#ifdef TEST_COVERAGE_MODE
112 long int k = 24288 * v_cl.getProcessingUnits();
113#else
114 long int k = 524288 * v_cl.getProcessingUnits();
115#endif
116
117 long int big_step = k / 4;
118 big_step = (big_step == 0)?1:big_step;
119
120 print_test_v( "Testing 2D vector with sfc curve reordering k<=",k);
121
122 // 2D test
123 for ( ; k >= 2 ; k-= decrement(k,big_step) )
124 {
125 BOOST_TEST_CHECKPOINT( "Testing 2D vector with sfc curve reordering k=" << k );
126
127 Box<2,float> box({0.0,0.0},{1.0,1.0});
128
129 // Boundary conditions
130 size_t bc[2]={NON_PERIODIC,NON_PERIODIC};
131
132 vector_dist<2,float, Point_test<float> > vd(k,box,bc,Ghost<2,float>(0.01));
133
134 auto it = vd.getIterator();
135
136 while (it.isNext())
137 {
138 auto key = it.get();
139
140 vd.getPos(key)[0] = ud(eg);
141 vd.getPos(key)[1] = ud(eg);
142
143 ++it;
144 }
145
146 vd.map();
147
148 // in case of SE_CLASS3 get a cell-list without ghost get is an error
149
150 // Create first cell list
151
152 auto NN1 = vd.getCellList(0.01,true);
153
154 //Reorder a vector
155 vd.reorder_rcut(0.01);
156
157 // Create second cell list
158 auto NN2 = vd.getCellList(0.01,true);
159
160 //Check equality of cell sizes
161 for (size_t i = 0 ; i < NN1.getGrid().size() ; i++)
162 {
163 size_t n1 = NN1.getNelements(i);
164 size_t n2 = NN2.getNelements(i);
165
166 BOOST_REQUIRE_EQUAL(n1,n2);
167 }
168 }
169}
170
171BOOST_AUTO_TEST_SUITE( vector_dist_cell_list_test_suite )
172
173BOOST_AUTO_TEST_CASE( vector_dist_reorder_2d_test )
174{
175 test_reorder_sfc(reorder_opt::HILBERT);
176 test_reorder_sfc(reorder_opt::LINEAR);
177}
178
179BOOST_AUTO_TEST_CASE( vector_dist_reorder_cl_test )
180{
181 test_reorder_cl();
182}
183
184BOOST_AUTO_TEST_CASE( vector_dist_cl_random_vs_hilb_forces_test )
185{
186 Vcluster<> & v_cl = create_vcluster();
187
188 if (v_cl.getProcessingUnits() > 48)
189 {return;}
190
191 ///////////////////// INPUT DATA //////////////////////
192
193 // Dimensionality of the space
194 const size_t dim = 3;
195 // Cut-off radiuses. Can be put different number of values
196 openfpm::vector<float> cl_r_cutoff {0.05};
197 // The starting amount of particles (remember that this number is multiplied by number of processors you use for testing)
198 size_t cl_k_start = 10000;
199 // The lower threshold for number of particles
200 size_t cl_k_min = 1000;
201 // Ghost part of distributed vector
202 double ghost_part = 0.05;
203
204 ///////////////////////////////////////////////////////
205
206 //For different r_cut
207 for (size_t r = 0; r < cl_r_cutoff.size(); r++ )
208 {
209 //Cut-off radius
210 float r_cut = cl_r_cutoff.get(r);
211
212 //Number of particles
213 size_t k = cl_k_start * v_cl.getProcessingUnits();
214
215 std::string str("Testing " + std::to_string(dim) + "D vector's forces (random vs hilb celllist) k<=");
216
217 print_test_v(str,k);
218
219 //For different number of particles
220 for (size_t k_int = k ; k_int >= cl_k_min ; k_int/=2 )
221 {
222 BOOST_TEST_CHECKPOINT( "Testing " << dim << "D vector's forces (random vs hilb celllist) k<=" << k_int );
223
224 Box<dim,float> box;
225
226 for (size_t i = 0; i < dim; i++)
227 {
228 box.setLow(i,0.0);
229 box.setHigh(i,1.0);
230 }
231
232 // Boundary conditions
233 size_t bc[dim];
234
235 for (size_t i = 0; i < dim; i++)
236 bc[i] = PERIODIC;
237
238 vector_dist<dim,float, aggregate<float[dim]> > vd(k_int,box,bc,Ghost<dim,float>(ghost_part));
239
240 vector_dist<dim,float, aggregate<float[dim]> > vd2(k_int,box,bc,Ghost<dim,float>(ghost_part));
241
242 // Initialize dist vectors
243 vd_initialize_double<dim>(vd, vd2, v_cl, k_int);
244
245 vd.ghost_get<0>();
246 vd2.ghost_get<0>();
247
248 //Get a cell list
249
250 auto NN = vd.getCellList(r_cut);
251
252 //Calculate forces
253
254 calc_forces<dim>(NN,vd,r_cut);
255
256 //Get a cell list hilb
257
258 auto NN_hilb = vd2.getCellList_hilb(r_cut);
259
260 //Calculate forces
261 calc_forces_hilb<dim>(NN_hilb,vd2,r_cut);
262
263 // Calculate average
264 size_t count = 1;
265 Point<dim,float> avg;
266 for (size_t i = 0 ; i < dim ; i++) {avg.get(i) = 0.0;}
267
268 auto it_v2 = vd.getIterator();
269 while (it_v2.isNext())
270 {
271 //key
272 vect_dist_key_dx key = it_v2.get();
273
274 for (size_t i = 0; i < dim; i++)
275 {avg.get(i) += fabs(vd.getProp<0>(key)[i]);}
276
277 ++count;
278 ++it_v2;
279 }
280
281 for (size_t i = 0 ; i < dim ; i++) {avg.get(i) /= count;}
282
283 auto it_v = vd.getIterator();
284 while (it_v.isNext())
285 {
286 //key
287 vect_dist_key_dx key = it_v.get();
288
289 for (size_t i = 0; i < dim; i++)
290 {
291 auto a1 = vd.getProp<0>(key)[i];
292 auto a2 = vd2.getProp<0>(key)[i];
293
294 //Check that the forces are (almost) equal
295 float per = 0.1;
296 if (a1 != 0.0)
297 per = fabs(0.1*avg.get(i)/a1);
298
299 BOOST_REQUIRE_CLOSE((float)a1,(float)a2,per);
300 }
301
302 ++it_v;
303 }
304 }
305 }
306}
307
308BOOST_AUTO_TEST_CASE( vector_dist_cl_random_vs_reorder_forces_test )
309{
310 Vcluster<> & v_cl = create_vcluster();
311
312 if (v_cl.getProcessingUnits() > 48)
313 return;
314
315 ///////////////////// INPUT DATA //////////////////////
316
317 // Dimensionality of the space
318 const size_t dim = 3;
319 // Cut-off radiuses. Can be put different number of values
320 openfpm::vector<float> cl_r_cutoff {0.01};
321 // The starting amount of particles (remember that this number is multiplied by number of processors you use for testing)
322 size_t cl_k_start = 10000;
323 // The lower threshold for number of particles
324 size_t cl_k_min = 1000;
325 // Ghost part of distributed vector
326 double ghost_part = 0.01;
327
328 ///////////////////////////////////////////////////////
329
330 //For different r_cut
331 for (size_t r = 0; r < cl_r_cutoff.size(); r++ )
332 {
333 //Cut-off radius
334 float r_cut = cl_r_cutoff.get(r);
335
336 //Number of particles
337 size_t k = cl_k_start * v_cl.getProcessingUnits();
338
339 std::string str("Testing " + std::to_string(dim) + "D vector's forces (random vs reorder) k<=");
340
341 print_test_v(str,k);
342
343 //For different number of particles
344 for (size_t k_int = k ; k_int >= cl_k_min ; k_int/=2 )
345 {
346 BOOST_TEST_CHECKPOINT( "Testing " << dim << "D vector's forces (random vs reorder) k<=" << k_int );
347
348 Box<dim,float> box;
349
350 for (size_t i = 0; i < dim; i++)
351 {
352 box.setLow(i,0.0);
353 box.setHigh(i,1.0);
354 }
355
356 // Boundary conditions
357 size_t bc[dim];
358
359 for (size_t i = 0; i < dim; i++)
360 bc[i] = PERIODIC;
361
362 vector_dist<dim,float, aggregate<float[dim], float[dim]> > vd(k_int,box,bc,Ghost<dim,float>(ghost_part));
363
364 // Initialize vd
365 vd_initialize<dim,decltype(vd)>(vd, v_cl);
366
367 vd.ghost_get<0>();
368
369 //Get a cell list
370
371 auto NN1 = vd.getCellList(r_cut);
372
373 //Calculate forces '0'
374
375 calc_forces<dim>(NN1,vd,r_cut);
376
377 //Reorder and get a cell list again
378
379 vd.reorder(4);
380
381 vd.ghost_get<0>();
382
383 auto NN2 = vd.getCellList(r_cut);
384
385 //Calculate forces '1'
386 calc_forces<dim,1>(NN2,vd,r_cut);
387
388 // Calculate average (For Coverty scan we start from 1)
389 size_t count = 1;
390 Point<dim,float> avg;
391 for (size_t i = 0 ; i < dim ; i++) {avg.get(i) = 0.0;}
392
393 auto it_v2 = vd.getIterator();
394 while (it_v2.isNext())
395 {
396 //key
397 vect_dist_key_dx key = it_v2.get();
398
399 for (size_t i = 0; i < dim; i++)
400 avg.get(i) += fabs(vd.getProp<0>(key)[i]);
401
402 ++count;
403 ++it_v2;
404 }
405
406 for (size_t i = 0 ; i < dim ; i++) {avg.get(i) /= count;}
407
408 //Test for equality of forces
409 auto it_v = vd.getDomainIterator();
410
411 while (it_v.isNext())
412 {
413 //key
414 vect_dist_key_dx key = it_v.get();
415
416 for (size_t i = 0; i < dim; i++)
417 {
418 float a1 = vd.getProp<0>(key)[i];
419 float a2 = vd.getProp<1>(key)[i];
420
421 //Check that the forces are (almost) equal
422 float per = 0.1;
423 if (a1 != 0.0)
424 per = fabs(0.1*avg.get(i)/a1);
425
426 BOOST_REQUIRE_CLOSE(a1,a2,per);
427 }
428
429 ++it_v;
430 }
431 }
432 }
433}
434
435BOOST_AUTO_TEST_CASE( vector_dist_symmetric_cell_list )
436{
437 Vcluster<> & v_cl = create_vcluster();
438
439 if (v_cl.getProcessingUnits() > 24)
440 return;
441
442 float L = 1000.0;
443
444 // set the seed
445 // create the random generator engine
446 std::srand(0);
447 std::default_random_engine eg;
448 std::uniform_real_distribution<float> ud(-L,L);
449
450 long int k = 4096 * v_cl.getProcessingUnits();
451
452 long int big_step = k / 4;
453 big_step = (big_step == 0)?1:big_step;
454
455 print_test_v("Testing 3D periodic vector symmetric cell-list k=",k);
456 BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list k=" << k );
457
458 Box<3,float> box({-L,-L,-L},{L,L,L});
459
460 // Boundary conditions
461 size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
462
463 float r_cut = 100.0;
464
465 // ghost
466 Ghost<3,float> ghost(r_cut);
467
468 // Point and global id
469 struct point_and_gid
470 {
471 size_t id;
472 Point<3,float> xq;
473
474 bool operator<(const struct point_and_gid & pag) const
475 {
476 return (id < pag.id);
477 }
478 };
479
480 typedef aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
481
482 // Distributed vector
483 vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_GHOST);
484 size_t start = vd.init_size_accum(k);
485
486 auto it = vd.getIterator();
487
488 while (it.isNext())
489 {
490 auto key = it.get();
491
492 vd.getPosWrite(key)[0] = ud(eg);
493 vd.getPosWrite(key)[1] = ud(eg);
494 vd.getPosWrite(key)[2] = ud(eg);
495
496 // Fill some properties randomly
497
498 vd.getPropWrite<0>(key) = 0;
499 vd.getPropWrite<1>(key) = 0;
500 vd.getPropWrite<2>(key) = key.getKey() + start;
501
502 ++it;
503 }
504
505 vd.map();
506
507 // sync the ghost
508 vd.ghost_get<0,2>();
509
510 auto NN = vd.getCellList(r_cut);
511 auto p_it = vd.getDomainIterator();
512
513 while (p_it.isNext())
514 {
515 auto p = p_it.get();
516
517 Point<3,float> xp = vd.getPosRead(p);
518
519 auto Np = NN.getNNIterator(NN.getCell(xp));
520
521 while (Np.isNext())
522 {
523 auto q = Np.get();
524
525 if (p.getKey() == q)
526 {
527 ++Np;
528 continue;
529 }
530
531 // repulsive
532
533 Point<3,float> xq = vd.getPosRead(q);
534 Point<3,float> f = (xp - xq);
535
536 float distance = f.norm();
537
538 // Particle should be inside 2 * r_cut range
539
540 if (distance < r_cut )
541 {
542 vd.getPropWrite<0>(p)++;
543 vd.getPropWrite<3>(p).add();
544 vd.getPropWrite<3>(p).last().xq = xq;
545 vd.getPropWrite<3>(p).last().id = vd.getPropWrite<2>(q);
546 }
547
548 ++Np;
549 }
550
551 ++p_it;
552 }
553
554 // We now try symmetric Cell-list
555
556 auto NN2 = vd.getCellListSym(r_cut);
557
558 auto p_it2 = vd.getDomainIterator();
559
560 while (p_it2.isNext())
561 {
562 auto p = p_it2.get();
563
564 Point<3,float> xp = vd.getPosRead(p);
565
566 auto Np = NN2.getNNIteratorSym<NO_CHECK>(NN2.getCell(xp),p.getKey(),vd.getPosVector());
567
568 while (Np.isNext())
569 {
570 auto q = Np.get();
571
572 if (p.getKey() == q)
573 {
574 ++Np;
575 continue;
576 }
577
578 // repulsive
579
580 Point<3,float> xq = vd.getPosRead(q);
581 Point<3,float> f = (xp - xq);
582
583 float distance = f.norm();
584
585 // Particle should be inside r_cut range
586
587 if (distance < r_cut )
588 {
589 vd.getPropWrite<1>(p)++;
590 vd.getPropWrite<1>(q)++;
591
592 vd.getPropWrite<4>(p).add();
593 vd.getPropWrite<4>(q).add();
594
595 vd.getPropWrite<4>(p).last().xq = xq;
596 vd.getPropWrite<4>(q).last().xq = xp;
597 vd.getPropWrite<4>(p).last().id = vd.getProp<2>(q);
598 vd.getPropWrite<4>(q).last().id = vd.getProp<2>(p);
599 }
600
601 ++Np;
602 }
603
604 ++p_it2;
605 }
606
607 vd.ghost_put<add_,1>();
608 vd.ghost_put<merge_,4>();
609
610 auto p_it3 = vd.getDomainIterator();
611
612 bool ret = true;
613 while (p_it3.isNext())
614 {
615 auto p = p_it3.get();
616
617 ret &= vd.getPropRead<1>(p) == vd.getPropRead<0>(p);
618
619 vd.getPropWrite<3>(p).sort();
620 vd.getPropWrite<4>(p).sort();
621
622 ret &= vd.getPropRead<3>(p).size() == vd.getPropRead<4>(p).size();
623
624 for (size_t i = 0 ; i < vd.getPropRead<3>(p).size() ; i++)
625 ret &= vd.getPropRead<3>(p).get(i).id == vd.getPropRead<4>(p).get(i).id;
626
627 if (ret == false)
628 {
629 std::cout << vd.getPropRead<3>(p).size() << " " << vd.getPropRead<4>(p).size() << std::endl;
630
631 for (size_t i = 0 ; i < vd.getPropRead<3>(p).size() ; i++)
632 std::cout << vd.getPropRead<3>(p).get(i).id << " " << vd.getPropRead<4>(p).get(i).id << std::endl;
633
634 std::cout << vd.getPropRead<1>(p) << " A " << vd.getPropRead<0>(p) << std::endl;
635
636 break;
637 }
638
639 ++p_it3;
640 }
641
642 BOOST_REQUIRE_EQUAL(ret,true);
643}
644
645BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_cell_list )
646{
647 Vcluster<> & v_cl = create_vcluster();
648
649 if (v_cl.getProcessingUnits() > 24)
650 return;
651
652 float L = 1000.0;
653
654 // set the seed
655 // create the random generator engine
656 std::default_random_engine eg(1132312*v_cl.getProcessUnitID());
657 std::uniform_real_distribution<float> ud(-L,L);
658
659 long int k = 4096 * v_cl.getProcessingUnits();
660
661 long int big_step = k / 4;
662 big_step = (big_step == 0)?1:big_step;
663
664 print_test_v("Testing 3D periodic vector symmetric crs cell-list k=",k);
665 BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric crs cell-list k=" << k );
666
667 Box<3,float> box({-L,-L,-L},{L,L,L});
668
669 // Boundary conditions
670 size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
671
672 float r_cut = 100.0;
673
674 // ghost
675 Ghost<3,float> ghost(r_cut);
676 Ghost<3,float> ghost2(r_cut);
677 ghost2.setLow(0,0.0);
678 ghost2.setLow(1,0.0);
679 ghost2.setLow(2,0.0);
680
681 // Point and global id
682 struct point_and_gid
683 {
684 size_t id;
685 Point<3,float> xq;
686
687 bool operator<(const struct point_and_gid & pag) const
688 {
689 return (id < pag.id);
690 }
691 };
692
693 typedef aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
694
695 // Distributed vector
696 vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_GHOST);
697 vector_dist<3,float, part_prop > vd2(k,box,bc,ghost2,BIND_DEC_TO_GHOST);
698 size_t start = vd.init_size_accum(k);
699
700 auto it = vd.getIterator();
701
702 while (it.isNext())
703 {
704 auto key = it.get();
705
706 vd.getPosWrite(key)[0] = ud(eg);
707 vd.getPosWrite(key)[1] = ud(eg);
708 vd.getPosWrite(key)[2] = ud(eg);
709
710 vd2.getPosWrite(key)[0] = vd.getPosRead(key)[0];
711 vd2.getPosWrite(key)[1] = vd.getPosRead(key)[1];
712 vd2.getPosWrite(key)[2] = vd.getPosRead(key)[2];
713
714 // Fill some properties randomly
715
716 vd.getPropWrite<0>(key) = 0;
717 vd.getPropWrite<1>(key) = 0;
718 vd.getPropWrite<2>(key) = key.getKey() + start;
719
720 vd2.getPropWrite<0>(key) = 0;
721 vd2.getPropWrite<1>(key) = 0;
722 vd2.getPropWrite<2>(key) = key.getKey() + start;
723
724 ++it;
725 }
726
727 vd.map();
728 vd2.map();
729
730 // sync the ghost
731 vd.ghost_get<0,2>();
732 vd2.ghost_get<0,2>();
733
734 auto NN = vd.getCellList(r_cut);
735 auto p_it = vd.getDomainIterator();
736
737 while (p_it.isNext())
738 {
739 auto p = p_it.get();
740
741 Point<3,float> xp = vd.getPosRead(p);
742
743 auto Np = NN.getNNIterator(NN.getCell(xp));
744
745 while (Np.isNext())
746 {
747 auto q = Np.get();
748
749 if (p.getKey() == q)
750 {
751 ++Np;
752 continue;
753 }
754
755 // repulsive
756
757 Point<3,float> xq = vd.getPosRead(q);
758 Point<3,float> f = (xp - xq);
759
760 float distance = f.norm();
761
762 // Particle should be inside 2 * r_cut range
763
764 if (distance < r_cut )
765 {
766 vd.getPropWrite<0>(p)++;
767 vd.getPropWrite<3>(p).add();
768 vd.getPropWrite<3>(p).last().xq = xq;
769 vd.getPropWrite<3>(p).last().id = vd.getPropRead<2>(q);
770 }
771
772 ++Np;
773 }
774
775 ++p_it;
776 }
777
778 // We now try symmetric Cell-list
779
780 auto NN2 = vd2.getCellListSym(r_cut);
781
782 // In case of CRS we have to iterate particles within some cells
783 // here we define whichone
784 auto p_it2 = vd2.getParticleIteratorCRS_Cell(NN2);
785
786 // For each particle
787 while (p_it2.isNext())
788 {
789 auto p = p_it2.get();
790
791 Point<3,float> xp = vd2.getPosRead(p);
792
793 auto Np = p_it2.getNNIteratorCSR(vd2.getPosVector());
794
795 while (Np.isNext())
796 {
797 auto q = Np.get();
798
799 if (p == q)
800 {
801 ++Np;
802 continue;
803 }
804
805 // repulsive
806
807 Point<3,float> xq = vd2.getPosRead(q);
808 Point<3,float> f = (xp - xq);
809
810 float distance = f.norm();
811
812 // Particle should be inside r_cut range
813
814 if (distance < r_cut )
815 {
816 vd2.getPropWrite<1>(p)++;
817 vd2.getPropWrite<1>(q)++;
818
819 vd2.getPropWrite<4>(p).add();
820 vd2.getPropWrite<4>(q).add();
821
822 vd2.getPropWrite<4>(p).last().xq = xq;
823 vd2.getPropWrite<4>(q).last().xq = xp;
824 vd2.getPropWrite<4>(p).last().id = vd2.getPropRead<2>(q);
825 vd2.getPropWrite<4>(q).last().id = vd2.getPropRead<2>(p);
826 }
827
828 ++Np;
829 }
830
831 ++p_it2;
832 }
833
834 vd2.ghost_put<add_,1>(NO_CHANGE_ELEMENTS);
835 vd2.ghost_put<merge_,4>();
836
837#ifdef SE_CLASS3
838 vd2.getDomainIterator();
839#endif
840
841 auto p_it3 = vd.getDomainIterator();
842
843 bool ret = true;
844 while (p_it3.isNext())
845 {
846 auto p = p_it3.get();
847
848 ret &= vd2.getPropRead<1>(p) == vd.getPropRead<0>(p);
849
850
851 vd.getPropWrite<3>(p).sort();
852 vd2.getPropWrite<4>(p).sort();
853
854 ret &= vd.getPropRead<3>(p).size() == vd2.getPropRead<4>(p).size();
855
856 for (size_t i = 0 ; i < vd.getPropRead<3>(p).size() ; i++)
857 ret &= vd.getPropRead<3>(p).get(i).id == vd2.getPropRead<4>(p).get(i).id;
858
859 if (ret == false)
860 break;
861
862 ++p_it3;
863 }
864
865 BOOST_REQUIRE_EQUAL(ret,true);
866}
867
868template<typename VerletList>
869void test_vd_symmetric_verlet_list()
870{
871 Vcluster<> & v_cl = create_vcluster();
872
873 if (v_cl.getProcessingUnits() > 24)
874 return;
875
876 float L = 1000.0;
877
878 // set the seed
879 // create the random generator engine
880 std::srand(0);
881 std::default_random_engine eg;
882 std::uniform_real_distribution<float> ud(-L,L);
883
884 long int k = 4096 * v_cl.getProcessingUnits();
885
886 long int big_step = k / 4;
887 big_step = (big_step == 0)?1:big_step;
888
889 print_test_v("Testing 3D periodic vector symmetric cell-list k=",k);
890 BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric verlet-list k=" << k );
891
892 Box<3,float> box({-L,-L,-L},{L,L,L});
893
894 // Boundary conditions
895 size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
896
897 float r_cut = 100.0;
898
899 // ghost
900 Ghost<3,float> ghost(r_cut);
901
902 // Point and global id
903 struct point_and_gid
904 {
905 size_t id;
906 Point<3,float> xq;
907
908 bool operator<(const struct point_and_gid & pag) const
909 {
910 return (id < pag.id);
911 }
912 };
913
914 typedef aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
915
916 // Distributed vector
917 vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_GHOST);
918 size_t start = vd.init_size_accum(k);
919
920 auto it = vd.getIterator();
921
922 while (it.isNext())
923 {
924 auto key = it.get();
925
926 vd.getPosWrite(key)[0] = ud(eg);
927 vd.getPosWrite(key)[1] = ud(eg);
928 vd.getPosWrite(key)[2] = ud(eg);
929
930 // Fill some properties randomly
931
932 vd.template getPropWrite<0>(key) = 0;
933 vd.template getPropWrite<1>(key) = 0;
934 vd.template getPropWrite<2>(key) = key.getKey() + start;
935
936 ++it;
937 }
938
939 vd.map();
940
941 // sync the ghost
942 vd.template ghost_get<0,2>();
943
944 auto NN = vd.template getVerlet<VerletList>(r_cut);
945 auto p_it = vd.getDomainIterator();
946
947 while (p_it.isNext())
948 {
949 auto p = p_it.get();
950
951 Point<3,float> xp = vd.getPosRead(p);
952
953 auto Np = NN.getNNIterator(p.getKey());
954
955 while (Np.isNext())
956 {
957 auto q = Np.get();
958
959 if (p.getKey() == q)
960 {
961 ++Np;
962 continue;
963 }
964
965 // repulsive
966
967 Point<3,float> xq = vd.getPosRead(q);
968 Point<3,float> f = (xp - xq);
969
970 float distance = f.norm();
971
972 // Particle should be inside 2 * r_cut range
973
974 if (distance < r_cut )
975 {
976 vd.template getPropWrite<0>(p)++;
977 vd.template getPropWrite<3>(p).add();
978 vd.template getPropWrite<3>(p).last().xq = xq;
979 vd.template getPropWrite<3>(p).last().id = vd.template getPropRead<2>(q);
980 }
981
982 ++Np;
983 }
984
985 ++p_it;
986 }
987
988 // We now try symmetric Cell-list
989
990 auto NN2 = vd.template getVerletSym<VerletList>(r_cut);
991
992 auto p_it2 = vd.getDomainIterator();
993
994 while (p_it2.isNext())
995 {
996 auto p = p_it2.get();
997
998 Point<3,float> xp = vd.getPosRead(p);
999
1000 auto Np = NN2.template getNNIterator<NO_CHECK>(p.getKey());
1001
1002 while (Np.isNext())
1003 {
1004 auto q = Np.get();
1005
1006 if (p.getKey() == q)
1007 {
1008 ++Np;
1009 continue;
1010 }
1011
1012 // repulsive
1013
1014 Point<3,float> xq = vd.getPosRead(q);
1015 Point<3,float> f = (xp - xq);
1016
1017 float distance = f.norm();
1018
1019 // Particle should be inside r_cut range
1020
1021 if (distance < r_cut )
1022 {
1023 vd.template getPropWrite<1>(p)++;
1024 vd.template getPropWrite<1>(q)++;
1025
1026 vd.template getPropWrite<4>(p).add();
1027 vd.template getPropWrite<4>(q).add();
1028
1029 vd.template getPropWrite<4>(p).last().xq = xq;
1030 vd.template getPropWrite<4>(q).last().xq = xp;
1031 vd.template getPropWrite<4>(p).last().id = vd.template getPropRead<2>(q);
1032 vd.template getPropWrite<4>(q).last().id = vd.template getPropRead<2>(p);
1033 }
1034
1035 ++Np;
1036 }
1037
1038 ++p_it2;
1039 }
1040
1041 vd.template ghost_put<add_,1>();
1042 vd.template ghost_put<merge_,4>();
1043
1044 auto p_it3 = vd.getDomainIterator();
1045
1046 bool ret = true;
1047 while (p_it3.isNext())
1048 {
1049 auto p = p_it3.get();
1050
1051 ret &= vd.template getPropRead<1>(p) == vd.template getPropRead<0>(p);
1052
1053 vd.template getPropWrite<3>(p).sort();
1054 vd.template getPropWrite<4>(p).sort();
1055
1056 ret &= vd.template getPropRead<3>(p).size() == vd.template getPropRead<4>(p).size();
1057
1058 for (size_t i = 0 ; i < vd.template getPropRead<3>(p).size() ; i++)
1059 ret &= vd.template getPropRead<3>(p).get(i).id == vd.template getPropRead<4>(p).get(i).id;
1060
1061 if (ret == false)
1062 break;
1063
1064 ++p_it3;
1065 }
1066
1067 BOOST_REQUIRE_EQUAL(ret,true);
1068}
1069
1070BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list )
1071{
1072 test_vd_symmetric_verlet_list<VERLET_MEMFAST(3,float)>();
1073 test_vd_symmetric_verlet_list<VERLET_MEMBAL(3,float)>();
1074 test_vd_symmetric_verlet_list<VERLET_MEMMW(3,float)>();
1075}
1076
1077template<typename VerletList>
1078void vector_sym_verlet_list_nb()
1079{
1080 Vcluster<> & v_cl = create_vcluster();
1081
1082 if (v_cl.getProcessingUnits() > 24)
1083 return;
1084
1085 float L = 1000.0;
1086
1087 // set the seed
1088 // create the random generator engine
1089 std::srand(0);
1090 std::default_random_engine eg;
1091 std::uniform_real_distribution<float> ud(-L,L);
1092
1093 long int k = 4096 * v_cl.getProcessingUnits();
1094
1095 long int big_step = k / 4;
1096 big_step = (big_step == 0)?1:big_step;
1097
1098 print_test_v("Testing 3D periodic vector symmetric cell-list no bottom k=",k);
1099 BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list no bottom k=" << k );
1100
1101 Box<3,float> box({-L,-L,-L},{L,L,L});
1102
1103 // Boundary conditions
1104 size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
1105
1106 float r_cut = 100.0;
1107
1108 // ghost
1109 Ghost<3,float> ghost(r_cut);
1110 Ghost<3,float> ghost2(r_cut);
1111 ghost2.setLow(2,0.0);
1112
1113 // Point and global id
1114 struct point_and_gid
1115 {
1116 size_t id;
1117 Point<3,float> xq;
1118
1119 bool operator<(const struct point_and_gid & pag) const
1120 {
1121 return (id < pag.id);
1122 }
1123 };
1124
1125 typedef aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
1126
1127 // 3D test
1128 for (size_t s = 0 ; s < 8 ; s++)
1129 {
1130
1131 // Distributed vector
1132 vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_GHOST);
1133 vector_dist<3,float, part_prop > vd2(k,box,bc,ghost2,BIND_DEC_TO_GHOST);
1134 size_t start = vd.init_size_accum(k);
1135
1136 auto it = vd.getIterator();
1137
1138 while (it.isNext())
1139 {
1140 auto key = it.get();
1141
1142 vd.getPosWrite(key)[0] = ud(eg);
1143 vd.getPosWrite(key)[1] = ud(eg);
1144 vd.getPosWrite(key)[2] = ud(eg);
1145
1146 vd2.getPosWrite(key)[0] = vd.getPosRead(key)[0];
1147 vd2.getPosWrite(key)[1] = vd.getPosRead(key)[1];
1148 vd2.getPosWrite(key)[2] = vd.getPosRead(key)[2];
1149
1150 // Fill some properties randomly
1151
1152 vd.template getPropWrite<0>(key) = 0;
1153 vd.template getPropWrite<1>(key) = 0;
1154 vd.template getPropWrite<2>(key) = key.getKey() + start;
1155
1156 vd2.template getPropWrite<0>(key) = 0;
1157 vd2.template getPropWrite<1>(key) = 0;
1158 vd2.template getPropWrite<2>(key) = key.getKey() + start;
1159
1160 ++it;
1161 }
1162
1163 vd.map();
1164 vd2.map();
1165
1166 // sync the ghost
1167 vd.template ghost_get<0,2>();
1168 vd2.template ghost_get<0,2>();
1169
1170 auto NN = vd.template getVerlet<VerletList>(r_cut);
1171 auto p_it = vd.getDomainIterator();
1172
1173 while (p_it.isNext())
1174 {
1175 auto p = p_it.get();
1176
1177 Point<3,float> xp = vd.getPosRead(p);
1178
1179 auto Np = NN.getNNIterator(p.getKey());
1180
1181 while (Np.isNext())
1182 {
1183 auto q = Np.get();
1184
1185 if (p.getKey() == q)
1186 {
1187 ++Np;
1188 continue;
1189 }
1190
1191 // repulsive
1192
1193 Point<3,float> xq = vd.getPosRead(q);
1194 Point<3,float> f = (xp - xq);
1195
1196 float distance = f.norm();
1197
1198 // Particle should be inside 2 * r_cut range
1199
1200 if (distance < r_cut )
1201 {
1202 vd.template getPropWrite<0>(p)++;
1203 vd.template getPropWrite<3>(p).add();
1204 vd.template getPropWrite<3>(p).last().xq = xq;
1205 vd.template getPropWrite<3>(p).last().id = vd.template getPropRead<2>(q);
1206 }
1207
1208 ++Np;
1209 }
1210
1211 ++p_it;
1212 }
1213
1214 // We now try symmetric Cell-list
1215
1216 auto NN2 = vd2.template getVerletSym<VerletList>(r_cut);
1217
1218 auto p_it2 = vd2.getDomainIterator();
1219
1220 while (p_it2.isNext())
1221 {
1222 auto p = p_it2.get();
1223
1224 Point<3,float> xp = vd2.getPosRead(p);
1225
1226 auto Np = NN2.template getNNIterator<NO_CHECK>(p.getKey());
1227
1228 while (Np.isNext())
1229 {
1230 auto q = Np.get();
1231
1232 if (p.getKey() == q)
1233 {
1234 ++Np;
1235 continue;
1236 }
1237
1238 // repulsive
1239
1240 Point<3,float> xq = vd2.getPosRead(q);
1241 Point<3,float> f = (xp - xq);
1242
1243 float distance = f.norm();
1244
1245 // Particle should be inside r_cut range
1246
1247 if (distance < r_cut )
1248 {
1249 vd2.template getPropWrite<1>(p)++;
1250 vd2.template getPropWrite<1>(q)++;
1251
1252 vd2.template getPropWrite<4>(p).add();
1253 vd2.template getPropWrite<4>(q).add();
1254
1255 vd2.template getPropWrite<4>(p).last().xq = xq;
1256 vd2.template getPropWrite<4>(q).last().xq = xp;
1257 vd2.template getPropWrite<4>(p).last().id = vd2.template getPropRead<2>(q);
1258 vd2.template getPropWrite<4>(q).last().id = vd2.template getPropRead<2>(p);
1259 }
1260
1261 ++Np;
1262 }
1263
1264
1265 ++p_it2;
1266 }
1267
1268 vd2.template ghost_put<add_,1>();
1269 vd2.template ghost_put<merge_,4>();
1270
1271#ifdef SE_CLASS3
1272 vd2.getDomainIterator();
1273#endif
1274
1275 auto p_it3 = vd.getDomainIterator();
1276
1277 bool ret = true;
1278 while (p_it3.isNext())
1279 {
1280 auto p = p_it3.get();
1281
1282 ret &= vd2.template getPropRead<1>(p) == vd.template getPropRead<0>(p);
1283
1284 vd.template getPropWrite<3>(p).sort();
1285 vd2.template getPropWrite<4>(p).sort();
1286
1287 ret &= vd.template getPropRead<3>(p).size() == vd2.template getPropRead<4>(p).size();
1288
1289 for (size_t i = 0 ; i < vd.template getPropRead<3>(p).size() ; i++)
1290 ret &= vd.template getPropRead<3>(p).get(i).id == vd2.template getPropRead<4>(p).get(i).id;
1291
1292 if (ret == false)
1293 break;
1294
1295 ++p_it3;
1296 }
1297
1298 BOOST_REQUIRE_EQUAL(ret,true);
1299 }
1300}
1301
1302BOOST_AUTO_TEST_CASE( vector_dist_symmetric_verlet_list_no_bottom )
1303{
1304 vector_sym_verlet_list_nb<VERLET_MEMFAST(3,float)>();
1305 vector_sym_verlet_list_nb<VERLET_MEMBAL(3,float)>();
1306 vector_sym_verlet_list_nb<VERLET_MEMMW(3,float)>();
1307
1308 vector_sym_verlet_list_nb<VERLET_MEMFAST_INT(3,float)>();
1309 vector_sym_verlet_list_nb<VERLET_MEMBAL_INT(3,float)>();
1310 vector_sym_verlet_list_nb<VERLET_MEMMW_INT(3,float)>();
1311}
1312
1313template<typename VerletList, typename part_prop> void test_crs_full(vector_dist<3,float, part_prop > & vd,
1314 vector_dist<3,float, part_prop > & vd2,
1315 std::default_random_engine & eg,
1316 std::uniform_real_distribution<float> & ud,
1317 size_t start,
1318 float r_cut)
1319{
1320 auto it = vd.getIterator();
1321
1322 while (it.isNext())
1323 {
1324 auto key = it.get();
1325
1326 vd.getPosWrite(key)[0] = ud(eg);
1327 vd.getPosWrite(key)[1] = ud(eg);
1328 vd.getPosWrite(key)[2] = ud(eg);
1329
1330 vd2.getPosWrite(key)[0] = vd.getPosRead(key)[0];
1331 vd2.getPosWrite(key)[1] = vd.getPosRead(key)[1];
1332 vd2.getPosWrite(key)[2] = vd.getPosRead(key)[2];
1333
1334 // Fill some properties randomly
1335
1336 vd.template getPropWrite<0>(key) = 0;
1337 vd.template getPropWrite<1>(key) = 0;
1338 vd.template getPropWrite<2>(key) = key.getKey() + start;
1339
1340 vd2.template getPropWrite<0>(key) = 0;
1341 vd2.template getPropWrite<1>(key) = 0;
1342 vd2.template getPropWrite<2>(key) = key.getKey() + start;
1343
1344 ++it;
1345 }
1346
1347 vd.map();
1348 vd2.map();
1349
1350 // sync the ghost
1351 vd.template ghost_get<0,2>();
1352 vd2.template ghost_get<0,2>();
1353
1354 auto NN = vd.template getVerlet<VerletList>(r_cut);
1355 auto p_it = vd.getDomainIterator();
1356
1357 while (p_it.isNext())
1358 {
1359 auto p = p_it.get();
1360
1361 Point<3,float> xp = vd.getPosRead(p);
1362
1363 auto Np = NN.getNNIterator(p.getKey());
1364
1365 while (Np.isNext())
1366 {
1367 auto q = Np.get();
1368
1369 if (p.getKey() == q)
1370 {
1371 ++Np;
1372 continue;
1373 }
1374
1375 // repulsive
1376
1377 Point<3,float> xq = vd.getPosRead(q);
1378 Point<3,float> f = (xp - xq);
1379
1380 float distance = f.norm();
1381
1382 // Particle should be inside 2 * r_cut range
1383
1384 if (distance < r_cut )
1385 {
1386 vd.template getPropWrite<0>(p)++;
1387 vd.template getPropWrite<3>(p).add();
1388 vd.template getPropWrite<3>(p).last().xq = xq;
1389 vd.template getPropWrite<3>(p).last().id = vd.template getPropRead<2>(q);
1390 }
1391
1392 ++Np;
1393 }
1394
1395 ++p_it;
1396 }
1397
1398 // We now try symmetric Verlet-list Crs scheme
1399
1400 auto NN2 = vd2.template getVerletCrs<VerletList>(r_cut);
1401
1402 // Because iterating across particles in the CSR scheme require a Cell-list
1403 auto p_it2 = vd2.getParticleIteratorCRS_Cell(NN2.getInternalCellList());
1404
1405 while (p_it2.isNext())
1406 {
1407 auto p = p_it2.get();
1408
1409 Point<3,float> xp = vd2.getPosRead(p);
1410
1411 auto Np = NN2.template getNNIterator<NO_CHECK>(p);
1412
1413 while (Np.isNext())
1414 {
1415 auto q = Np.get();
1416
1417 if (p == q)
1418 {
1419 ++Np;
1420 continue;
1421 }
1422
1423 // repulsive
1424
1425 Point<3,float> xq = vd2.getPosRead(q);
1426 Point<3,float> f = (xp - xq);
1427
1428 float distance = f.norm();
1429
1430 if (distance < r_cut )
1431 {
1432 vd2.template getPropWrite<1>(p)++;
1433 vd2.template getPropWrite<1>(q)++;
1434
1435 vd2.template getPropWrite<4>(p).add();
1436 vd2.template getPropWrite<4>(q).add();
1437
1438 vd2.template getPropWrite<4>(p).last().xq = xq;
1439 vd2.template getPropWrite<4>(q).last().xq = xp;
1440 vd2.template getPropWrite<4>(p).last().id = vd2.template getPropRead<2>(q);
1441 vd2.template getPropWrite<4>(q).last().id = vd2.template getPropRead<2>(p);
1442 }
1443
1444 ++Np;
1445 }
1446
1447 ++p_it2;
1448 }
1449
1450 vd2.template ghost_put<add_,1>();
1451 vd2.template ghost_put<merge_,4>();
1452
1453#ifdef SE_CLASS3
1454 vd2.getDomainIterator();
1455#endif
1456
1457 auto p_it3 = vd.getDomainIterator();
1458
1459 bool ret = true;
1460 while (p_it3.isNext())
1461 {
1462 auto p = p_it3.get();
1463
1464 ret &= vd2.template getPropRead<1>(p) == vd.template getPropRead<0>(p);
1465
1466 if (ret == false)
1467 {
1468 Point<3,float> xp = vd2.getPosRead(p);
1469 std::cout << "ERROR " << vd2.template getPropWrite<1>(p) << " " << vd.template getPropWrite<0>(p) << " " << xp.toString() << std::endl;
1470 }
1471
1472 vd.template getPropWrite<3>(p).sort();
1473 vd2.template getPropWrite<4>(p).sort();
1474
1475 ret &= vd.template getPropRead<3>(p).size() == vd2.template getPropRead<4>(p).size();
1476
1477 for (size_t i = 0 ; i < vd.template getPropRead<3>(p).size() ; i++)
1478 ret &= vd.template getPropRead<3>(p).get(i).id == vd2.template getPropRead<4>(p).get(i).id;
1479
1480 if (ret == false)
1481 break;
1482
1483 ++p_it3;
1484 }
1485
1486 BOOST_REQUIRE_EQUAL(ret,true);
1487}
1488
1489template<typename VerletList>
1490void test_csr_verlet_list()
1491{
1492 Vcluster<> & v_cl = create_vcluster();
1493
1494 if (v_cl.getProcessingUnits() > 24)
1495 return;
1496
1497 float L = 1000.0;
1498
1499 // set the seed
1500 // create the random generator engine
1501 std::srand(0);
1502 std::default_random_engine eg;
1503 std::uniform_real_distribution<float> ud(-L,L);
1504
1505 long int k = 4096 * v_cl.getProcessingUnits();
1506
1507 long int big_step = k / 4;
1508 big_step = (big_step == 0)?1:big_step;
1509
1510 print_test_v("Testing 3D periodic vector symmetric cell-list k=",k);
1511 BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list k=" << k );
1512
1513 Box<3,float> box({-L,-L,-L},{L,L,L});
1514
1515 // Boundary conditions
1516 size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
1517
1518 float r_cut = 100.0;
1519
1520 // ghost
1521 Ghost<3,float> ghost(r_cut);
1522 Ghost<3,float> ghost2(r_cut);
1523 ghost2.setLow(0,0.0);
1524 ghost2.setLow(1,0.0);
1525 ghost2.setLow(2,0.0);
1526
1527 // Point and global id
1528 struct point_and_gid
1529 {
1530 size_t id;
1531 Point<3,float> xq;
1532
1533 bool operator<(const struct point_and_gid & pag) const
1534 {
1535 return (id < pag.id);
1536 }
1537 };
1538
1539 typedef aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
1540
1541 // Distributed vector
1542 vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_GHOST);
1543 vector_dist<3,float, part_prop > vd2(k,box,bc,ghost2,BIND_DEC_TO_GHOST);
1544 size_t start = vd.init_size_accum(k);
1545
1546 test_crs_full<VerletList>(vd,vd2,eg,ud,start,r_cut);
1547}
1548
1549template<typename VerletList>
1550void test_csr_verlet_list_override()
1551{
1552 Vcluster<> & v_cl = create_vcluster();
1553
1554 if (v_cl.getProcessingUnits() > 24)
1555 return;
1556
1557 float L = 1000.0;
1558
1559 // set the seed
1560 // create the random generator engine
1561 std::srand(0);
1562 std::default_random_engine eg;
1563 std::uniform_real_distribution<float> ud(-L,L);
1564
1565 long int k = 4096 * v_cl.getProcessingUnits();
1566
1567 long int big_step = k / 4;
1568 big_step = (big_step == 0)?1:big_step;
1569
1570 print_test_v("Testing 3D periodic vector symmetric cell-list k=",k);
1571 BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list k=" << k );
1572
1573 Box<3,float> box({-L,-L,-L},{L,L,L});
1574
1575 // Boundary conditions
1576 size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
1577
1578 float r_cut = 100.0;
1579
1580 // ghost
1581 Ghost<3,float> ghost(r_cut);
1582 Ghost<3,float> ghost2(r_cut);
1583 ghost2.setLow(0,0.0);
1584 ghost2.setLow(1,0.0);
1585 ghost2.setLow(2,0.0);
1586
1587 // Point and global id
1588 struct point_and_gid
1589 {
1590 size_t id;
1591 Point<3,float> xq;
1592
1593 bool operator<(const struct point_and_gid & pag) const
1594 {
1595 return (id < pag.id);
1596 }
1597 };
1598
1599 typedef aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
1600
1601 size_t gdist_d[3];
1602 size_t gdist2_d[3];
1603
1604 gdist_d[0] = 1;
1605 gdist_d[1] = 2;
1606 gdist_d[2] = 5;
1607
1608 gdist2_d[0] = 1;
1609 gdist2_d[1] = 2;
1610 gdist2_d[2] = 5;
1611
1612 grid_sm<3,void> gdist(gdist_d);
1613 grid_sm<3,void> gdist2(gdist2_d);
1614
1615 // Distributed vector
1616 vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_GHOST,gdist_d);
1617 vector_dist<3,float, part_prop > vd2(k,box,bc,ghost2,BIND_DEC_TO_GHOST,gdist2_d);
1618 size_t start = vd.init_size_accum(k);
1619
1620 test_crs_full<VerletList>(vd,vd2,eg,ud,start,r_cut);
1621}
1622
1623BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list )
1624{
1625 test_csr_verlet_list<VERLET_MEMFAST(3,float)>();
1626 test_csr_verlet_list<VERLET_MEMBAL(3,float)>();
1627 test_csr_verlet_list<VERLET_MEMMW(3,float)>();
1628}
1629
1630BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list_dec_override )
1631{
1632 test_csr_verlet_list_override<VERLET_MEMFAST(3,float)>();
1633 test_csr_verlet_list_override<VERLET_MEMBAL(3,float)>();
1634 test_csr_verlet_list_override<VERLET_MEMMW(3,float)>();
1635}
1636
1637template <typename VerletList>
1638void test_vd_symmetric_crs_verlet()
1639{
1640 Vcluster<> & v_cl = create_vcluster();
1641
1642 if (v_cl.getProcessingUnits() > 24)
1643 return;
1644
1645 float L = 1000.0;
1646
1647 bool ret = true;
1648
1649 // set the seed
1650 // create the random generator engine
1651 std::srand(0);
1652 std::default_random_engine eg;
1653 std::uniform_real_distribution<float> ud(-L,L);
1654
1655 long int k = 4096 * v_cl.getProcessingUnits();
1656
1657 long int big_step = k / 4;
1658 big_step = (big_step == 0)?1:big_step;
1659
1660 print_test_v("Testing 3D periodic vector symmetric cell-list k=",k);
1661 BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list k=" << k );
1662
1663 Box<3,float> box({-L,-L,-L},{L,L,L});
1664
1665 // Boundary conditions
1666 size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
1667
1668 float r_cut = 100.0;
1669
1670 // ghost
1671 Ghost<3,float> ghost(r_cut);
1672 Ghost<3,float> ghost2(r_cut);
1673 ghost2.setLow(0,0.0);
1674 ghost2.setLow(1,0.0);
1675 ghost2.setLow(2,0.0);
1676
1677
1678 typedef aggregate<size_t> part_prop;
1679
1680 // Distributed vector
1681 vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_GHOST);
1682
1683 auto it = vd.getIterator();
1684
1685 while (it.isNext())
1686 {
1687 auto key = it.get();
1688
1689 vd.getPos(key)[0] = ud(eg);
1690 vd.getPos(key)[1] = ud(eg);
1691 vd.getPos(key)[2] = ud(eg);
1692
1693 // Fill some properties randomly
1694
1695 vd.getProp<0>(key) = 0;
1696
1697 ++it;
1698 }
1699
1700 vd.map();
1701
1702 // sync the ghost
1703 vd.ghost_get<0>();
1704
1705 // We now try symmetric Verlet-list Crs scheme
1706
1707 auto NN2 = vd.template getVerletCrs<VerletList>(r_cut);
1708
1709 // Because iterating across particles in the CSR scheme require a Cell-list
1710 auto p_it2 = vd.getParticleIteratorCRS_Cell(NN2.getInternalCellList());
1711 auto p_it3 = vd.getParticleIteratorCRS(NN2);
1712
1713 while (p_it2.isNext())
1714 {
1715 auto p = p_it2.get();
1716 auto p2 = p_it3.get();
1717
1718 ret &= (p == p2);
1719
1720 if (ret == false)
1721 break;
1722
1723 ++p_it2;
1724 ++p_it3;
1725 }
1726
1727 BOOST_REQUIRE_EQUAL(ret,true);
1728}
1729
1730BOOST_AUTO_TEST_CASE( vector_dist_symmetric_crs_verlet_list_partit )
1731{
1732 test_vd_symmetric_crs_verlet<VERLET_MEMFAST(3,float)>();
1733 test_vd_symmetric_crs_verlet<VERLET_MEMBAL(3,float)>();
1734 test_vd_symmetric_crs_verlet<VERLET_MEMMW(3,float)>();
1735}
1736
1737BOOST_AUTO_TEST_CASE( vector_dist_checking_unloaded_processors )
1738{
1739 Vcluster<> & v_cl = create_vcluster();
1740
1741 if (v_cl.getProcessingUnits() > 24)
1742 return;
1743
1744 float L = 200.0;
1745
1746 // set the seed
1747 // create the random generator engine
1748 std::srand(0);
1749 std::default_random_engine eg;
1750 std::uniform_real_distribution<float> ud(0,L);
1751
1752 long int k = 4096 * v_cl.getProcessingUnits();
1753
1754 long int big_step = k / 4;
1755 big_step = (big_step == 0)?1:big_step;
1756
1757 print_test_v("Testing 3D periodic vector symmetric cell-list (unload processors) k=",k);
1758 BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list (unload processors) k=" << k );
1759
1760 Box<3,float> box({0,0,0},{L,L,L});
1761
1762 // Boundary conditions
1763 size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
1764
1765 float r_cut = 100.0;
1766
1767 // ghost
1768 Ghost<3,float> ghost(r_cut);
1769 Ghost<3,float> ghost2(r_cut);
1770 ghost2.setLow(0,0.0);
1771 ghost2.setLow(1,0.0);
1772 ghost2.setLow(2,0.0);
1773
1774
1775 typedef aggregate<size_t> part_prop;
1776
1777 // Distributed vector
1778 vector_dist<3,float, part_prop > vd(k,box,bc,ghost,BIND_DEC_TO_GHOST);
1779
1780 auto it = vd.getIterator();
1781
1782 while (it.isNext())
1783 {
1784 auto key = it.get();
1785
1786 vd.getPos(key)[0] = ud(eg);
1787 vd.getPos(key)[1] = ud(eg);
1788 vd.getPos(key)[2] = ud(eg);
1789
1790 // Fill some properties randomly
1791
1792 vd.getProp<0>(key) = 0;
1793
1794 ++it;
1795 }
1796
1797 vd.map();
1798
1799 //
1800 if (v_cl.getProcessingUnits() >= 9)
1801 {
1802 size_t min = vd.size_local();
1803
1804 v_cl.min(min);
1805 v_cl.execute();
1806
1807 BOOST_REQUIRE_EQUAL(min,0ul);
1808 }
1809
1810
1811 // sync the ghost
1812 vd.ghost_get<0>();
1813
1814 //
1815 if (v_cl.getProcessingUnits() >= 9)
1816 {
1817 size_t min = vd.size_local_with_ghost() - vd.size_local();
1818
1819 v_cl.min(min);
1820 v_cl.execute();
1821
1822 BOOST_REQUIRE_EQUAL(min,0ul);
1823 }
1824}
1825
1826BOOST_AUTO_TEST_CASE( vector_dist_cell_list_multi_type )
1827{
1828 Vcluster<> & v_cl = create_vcluster();
1829
1830 if (v_cl.getProcessingUnits() > 24)
1831 return;
1832
1833 float L = 1000.0;
1834
1835 // set the seed
1836 // create the random generator engine
1837 std::srand(0);
1838 std::default_random_engine eg;
1839 std::uniform_real_distribution<float> ud(-L,L);
1840
1841 long int k = 4096 * v_cl.getProcessingUnits();
1842
1843 long int big_step = k / 4;
1844 big_step = (big_step == 0)?1:big_step;
1845
1846 print_test_v("Testing 3D periodic vector symmetric cell-list k=",k);
1847 BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list k=" << k );
1848
1849 Box<3,float> box({-L,-L,-L},{L,L,L});
1850
1851 // Boundary conditions
1852 size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
1853
1854 float r_cut = 100.0;
1855
1856 // ghost
1857 Ghost<3,float> ghost(r_cut);
1858
1859 typedef aggregate<size_t> part_prop;
1860
1861 // Distributed vector
1862 vector_dist<3,float, part_prop > vd(k,box,bc,ghost);
1863
1864 auto it = vd.getIterator();
1865
1866 while (it.isNext())
1867 {
1868 auto key = it.get();
1869
1870 vd.getPos(key)[0] = ud(eg);
1871 vd.getPos(key)[1] = ud(eg);
1872 vd.getPos(key)[2] = ud(eg);
1873
1874 ++it;
1875 }
1876
1877 vd.map();
1878
1879 // sync the ghost
1880 vd.ghost_get<0>();
1881
1882
1883 bool ret = true;
1884
1885 // We take different type of Cell-list
1886 auto NN = vd.getCellList<CELL_MEMFAST(3,float)>(r_cut);
1887 auto NN2 = vd.getCellList<CELL_MEMBAL(3,float)>(r_cut);
1888 auto NN3 = vd.getCellList<CELL_MEMMW(3,float)>(r_cut);
1889
1890 auto p_it = vd.getDomainIterator();
1891
1892 while (p_it.isNext())
1893 {
1894 auto p = p_it.get();
1895
1896 Point<3,float> xp = vd.getPos(p);
1897
1898 auto Np = NN.getNNIterator(NN.getCell(xp));
1899 auto Np2 = NN2.getNNIterator(NN2.getCell(xp));
1900 auto Np3 = NN3.getNNIterator(NN3.getCell(xp));
1901
1902 while (Np.isNext())
1903 {
1904 // first all cell-list must agree
1905
1906 ret &= (Np.isNext() == Np2.isNext()) && (Np3.isNext() == Np.isNext());
1907
1908 if (ret == false)
1909 break;
1910
1911 auto q = Np.get();
1912 auto q2 = Np2.get();
1913 auto q3 = Np3.get();
1914
1915 ret &= (q == q2) && (q == q3);
1916
1917 if (ret == false)
1918 break;
1919
1920 ++Np;
1921 ++Np2;
1922 ++Np3;
1923 }
1924
1925 ret &= (Np.isNext() == Np2.isNext()) && (Np.isNext() == Np3.isNext());
1926
1927 if (ret == false)
1928 break;
1929
1930 ++p_it;
1931 }
1932
1933 BOOST_REQUIRE_EQUAL(ret,true);
1934}
1935
1936// Point and global id
1937struct point_and_gid
1938{
1939 size_t id;
1940 Point<3,float> xq;
1941
1942 bool operator<(const struct point_and_gid & pag) const
1943 {
1944 return (id < pag.id);
1945 }
1946};
1947
1948template<typename vector_dist_mp>
1949void test_vector_dist_particle_NN_MP_iteration()
1950{
1951 Vcluster<> & v_cl = create_vcluster();
1952
1953 if (v_cl.getProcessingUnits() > 24)
1954 {return;}
1955
1956 float L = 1000.0;
1957
1958 // set the seed
1959 // create the random generator engine
1960 std::default_random_engine eg;
1961 eg.seed(v_cl.rank()*4533);
1962 std::uniform_real_distribution<float> ud(-L,L);
1963
1964 long int k = 4096 * v_cl.getProcessingUnits();
1965
1966 long int big_step = k / 4;
1967 big_step = (big_step == 0)?1:big_step;
1968
1969 print_test_v("Testing 3D periodic vector symmetric cell-list k=",k);
1970 BOOST_TEST_CHECKPOINT( "Testing 3D periodic vector symmetric cell-list k=" << k );
1971
1972 Box<3,float> box({-L,-L,-L},{L,L,L});
1973
1974 // Boundary conditions
1975 size_t bc[3]={PERIODIC,PERIODIC,PERIODIC};
1976
1977 float r_cut = 100.0;
1978
1979 // ghost
1980 Ghost<3,float> ghost(r_cut);
1981
1982// typedef aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
1983
1984 // Distributed vector
1985 vector_dist_mp vd(k,box,bc,ghost,BIND_DEC_TO_GHOST);
1986 size_t start = vd.init_size_accum(k);
1987
1988 auto it = vd.getIterator();
1989
1990 while (it.isNext())
1991 {
1992 auto key = it.get();
1993
1994 vd.getPosWrite(key)[0] = ud(eg);
1995 vd.getPosWrite(key)[1] = ud(eg);
1996 vd.getPosWrite(key)[2] = ud(eg);
1997
1998 // Fill some properties randomly
1999
2000 vd.template getPropWrite<0>(key) = 0;
2001 vd.template getPropWrite<1>(key) = 0;
2002 vd.template getPropWrite<2>(key) = key.getKey() + start;
2003
2004 ++it;
2005 }
2006
2007 vd.map();
2008
2009 // sync the ghost
2010 vd.template ghost_get<0,2>();
2011
2012 auto NN = vd.getCellList(r_cut);
2013 auto p_it = vd.getDomainIterator();
2014
2015 while (p_it.isNext())
2016 {
2017 auto p = p_it.get();
2018
2019 Point<3,float> xp = vd.getPosRead(p);
2020
2021 auto Np = NN.getNNIterator(NN.getCell(xp));
2022
2023 while (Np.isNext())
2024 {
2025 auto q = Np.get();
2026
2027 if (p.getKey() == q)
2028 {
2029 ++Np;
2030 continue;
2031 }
2032
2033 // repulsive
2034
2035 Point<3,float> xq = vd.getPosRead(q);
2036 Point<3,float> f = (xp - xq);
2037
2038 float distance = f.norm();
2039
2040 // Particle should be inside 2 * r_cut range
2041
2042 if (distance < r_cut )
2043 {
2044 vd.template getPropWrite<0>(p)++;
2045 vd.template getPropWrite<3>(p).add();
2046 vd.template getPropWrite<3>(p).last().xq = xq;
2047 vd.template getPropWrite<3>(p).last().id = vd.template getPropWrite<2>(q);
2048 }
2049
2050 ++Np;
2051 }
2052
2053 ++p_it;
2054 }
2055
2056 // We now divide the particles on 4 phases
2057
2058 openfpm::vector<vector_dist_mp> phases;
2059 phases.add( vector_dist_mp(vd.getDecomposition(),0));
2060 phases.add( vector_dist_mp(phases.get(0).getDecomposition(),0));
2061 phases.add( vector_dist_mp(phases.get(0).getDecomposition(),0));
2062 phases.add( vector_dist_mp(phases.get(0).getDecomposition(),0));
2063
2064 auto it2 = vd.getDomainIterator();
2065
2066 while (it2.isNext())
2067 {
2068 auto p = it2.get();
2069
2070 if (p.getKey() % 4 == 0)
2071 {
2072 phases.get(0).add();
2073 phases.get(0).getLastPos()[0] = vd.getPos(p)[0];
2074 phases.get(0).getLastPos()[1] = vd.getPos(p)[1];
2075 phases.get(0).getLastPos()[2] = vd.getPos(p)[2];
2076
2077 phases.get(0).template getLastProp<1>() = 0;
2078
2079 phases.get(0).template getLastProp<2>() = vd.template getProp<2>(p);
2080 }
2081 else if (p.getKey() % 4 == 1)
2082 {
2083 phases.get(1).add();
2084 phases.get(1).getLastPos()[0] = vd.getPos(p)[0];
2085 phases.get(1).getLastPos()[1] = vd.getPos(p)[1];
2086 phases.get(1).getLastPos()[2] = vd.getPos(p)[2];
2087
2088 phases.get(1).template getLastProp<1>() = 0;
2089
2090 phases.get(1).template getLastProp<2>() = vd.template getProp<2>(p);
2091 }
2092 else if (p.getKey() % 4 == 2)
2093 {
2094 phases.get(2).add();
2095 phases.get(2).getLastPos()[0] = vd.getPos(p)[0];
2096 phases.get(2).getLastPos()[1] = vd.getPos(p)[1];
2097 phases.get(2).getLastPos()[2] = vd.getPos(p)[2];
2098
2099 phases.get(2).template getLastProp<1>() = 0;
2100
2101 phases.get(2).template getLastProp<2>() = vd.template getProp<2>(p);
2102 }
2103 else
2104 {
2105 phases.get(3).add();
2106 phases.get(3).getLastPos()[0] = vd.getPos(p)[0];
2107 phases.get(3).getLastPos()[1] = vd.getPos(p)[1];
2108 phases.get(3).getLastPos()[2] = vd.getPos(p)[2];
2109
2110 phases.get(3).template getLastProp<1>() = 0;
2111
2112 phases.get(3).template getLastProp<2>() = vd.template getProp<2>(p);
2113 }
2114
2115 ++it2;
2116 }
2117
2118 // now we get all the Cell-lists
2119
2120 for (size_t i = 0 ; i < phases.size() ; i++)
2121 {
2122 phases.get(i).template ghost_get<0,1,2>();
2123 }
2124
2125 openfpm::vector<CellList<3, float, Mem_fast<>, shift<3, float> >> NN_ptr;
2126
2127 for (size_t i = 0 ; i < phases.size() ; i++)
2128 {
2129 NN_ptr.add(phases.get(i).getCellListSym(r_cut));
2130 }
2131
2132 // We now interact all the phases
2133
2134 for (size_t i = 0 ; i < phases.size() ; i++)
2135 {
2136 for (size_t j = 0 ; j < phases.size() ; j++)
2137 {
2138 auto p_it2 = phases.get(i).getDomainIterator();
2139
2140 while (p_it2.isNext())
2141 {
2142 auto p = p_it2.get();
2143
2144 Point<3,float> xp = phases.get(i).getPosRead(p);
2145
2146 auto Np = NN_ptr.get(j).getNNIteratorSymMP<NO_CHECK>(NN_ptr.get(j).getCell(xp),p.getKey(),phases.get(i).getPosVector(),phases.get(j).getPosVector());
2147
2148 while (Np.isNext())
2149 {
2150 auto q = Np.get();
2151
2152 if (p.getKey() == q && i == j)
2153 {
2154 ++Np;
2155 continue;
2156 }
2157
2158 // repulsive
2159
2160 Point<3,float> xq = phases.get(j).getPosRead(q);
2161 Point<3,float> f = (xp - xq);
2162
2163 float distance = f.norm();
2164
2165 // Particle should be inside r_cut range
2166
2167 if (distance < r_cut )
2168 {
2169 phases.get(i).template getPropWrite<1>(p)++;
2170 phases.get(j).template getPropWrite<1>(q)++;
2171
2172 phases.get(i).template getPropWrite<4>(p).add();
2173 phases.get(j).template getPropWrite<4>(q).add();
2174
2175 phases.get(i).template getPropWrite<4>(p).last().xq = xq;
2176 phases.get(j).template getPropWrite<4>(q).last().xq = xp;
2177 phases.get(i).template getPropWrite<4>(p).last().id = phases.get(j).template getProp<2>(q);
2178 phases.get(j).template getPropWrite<4>(q).last().id = phases.get(i).template getProp<2>(p);
2179 }
2180
2181 ++Np;
2182 }
2183
2184 ++p_it2;
2185 }
2186 }
2187 }
2188
2189 for (size_t i = 0 ; i < phases.size() ; i++)
2190 {
2191 phases.get(i).template ghost_put<add_,1>();
2192 phases.get(i).template ghost_put<merge_,4>();
2193 }
2194
2195 auto p_it3 = vd.getDomainIterator();
2196
2197 bool ret = true;
2198 while (p_it3.isNext())
2199 {
2200 auto p = p_it3.get();
2201
2202 int ph;
2203
2204 if (p.getKey() % 4 == 0)
2205 {ph = 0;}
2206 else if (p.getKey() % 4 == 1)
2207 {ph = 1;}
2208 else if (p.getKey() % 4 == 2)
2209 {ph = 2;}
2210 else
2211 {ph = 3;}
2212
2213 size_t pah = p.getKey()/4;
2214 ret &= phases.get(ph).template getPropRead<1>(pah) == vd.template getPropRead<0>(p);
2215
2216 vd.template getPropWrite<3>(p).sort();
2217 phases.get(ph).template getPropWrite<4>(pah).sort();
2218
2219 ret &= vd.template getPropRead<3>(p).size() == phases.get(ph).template getPropRead<4>(pah).size();
2220
2221 for (size_t i = 0 ; i < vd.template getPropRead<3>(p).size() ; i++)
2222 ret &= vd.template getPropRead<3>(p).get(i).id == phases.get(ph).template getPropRead<4>(pah).get(i).id;
2223
2224 if (ret == false)
2225 {
2226 std::cout << "Error on particle: " << vd.template getPropRead<2>(p) << " " << v_cl.rank() << std::endl;
2227
2228 std::cout << vd.template getPropRead<3>(p).size() << " " << phases.get(ph).template getPropRead<4>(pah).size() << " " << v_cl.rank() << std::endl;
2229
2230 for (size_t i = 0 ; i < vd.template getPropRead<3>(p).size() ; i++)
2231 std::cout << vd.template getPropRead<3>(p).get(i).id << " " << phases.get(ph).template getPropRead<4>(pah).get(i).id << " " << v_cl.rank() << std::endl;
2232
2233 std::cout << phases.get(ph).template getPropRead<1>(pah) << " A " << vd.template getPropRead<0>(p) << std::endl;
2234
2235 break;
2236 }
2237
2238 ++p_it3;
2239 }
2240
2241 BOOST_REQUIRE_EQUAL(ret,true);
2242}
2243
2244BOOST_AUTO_TEST_CASE( vector_dist_particle_NN_MP_iteration )
2245{
2246 typedef aggregate<size_t,size_t,size_t,openfpm::vector<point_and_gid>,openfpm::vector<point_and_gid>> part_prop;
2247
2248 test_vector_dist_particle_NN_MP_iteration<vector_dist<3,float, part_prop >>();
2249}
2250
2251BOOST_AUTO_TEST_SUITE_END()
2252