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 | |
15 | extern void print_test_v(std::string test, size_t sz); |
16 | extern long int decrement(long int k, long int step); |
17 | |
18 | ///////////////////////// test hilb /////////////////////////////// |
19 | |
20 | void 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 | |
98 | void 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 | |
171 | BOOST_AUTO_TEST_SUITE( vector_dist_cell_list_test_suite ) |
172 | |
173 | BOOST_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 | |
179 | BOOST_AUTO_TEST_CASE( vector_dist_reorder_cl_test ) |
180 | { |
181 | test_reorder_cl(); |
182 | } |
183 | |
184 | BOOST_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 | |
308 | BOOST_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 | |
435 | BOOST_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 | |
645 | BOOST_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 | |
868 | template<typename VerletList> |
869 | void 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 | |
1070 | BOOST_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 | |
1077 | template<typename VerletList> |
1078 | void 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 | |
1302 | BOOST_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 | |
1313 | template<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 | |
1489 | template<typename VerletList> |
1490 | void 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 | |
1549 | template<typename VerletList> |
1550 | void 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 | |
1623 | BOOST_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 | |
1630 | BOOST_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 | |
1637 | template <typename VerletList> |
1638 | void 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 | |
1730 | BOOST_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 | |
1737 | BOOST_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 | |
1826 | BOOST_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 |
1937 | struct 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 | |
1948 | template<typename vector_dist_mp> |
1949 | void 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 | |
2244 | BOOST_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 | |
2251 | BOOST_AUTO_TEST_SUITE_END() |
2252 | |