1/*
2 * SparseGrid_unit_tests.cpp
3 *
4 * Created on: Oct 22, 2017
5 * Author: i-bird
6 */
7
8#define DISABLE_MPI_WRITTERS
9
10#define BOOST_TEST_DYN_LINK
11#include <boost/test/unit_test.hpp>
12#include "SparseGrid/SparseGrid.hpp"
13#include "NN/CellList/CellDecomposer.hpp"
14#include <math.h>
15//#include "util/debug.hpp"
16
17BOOST_AUTO_TEST_SUITE( sparse_grid_test )
18
19template <typename grid_type, typename cell_decomposer>
20size_t fill_sphere(grid_type & grid, cell_decomposer & cdsm)
21{
22 size_t tot_count = 0;
23 double r = 0.3;
24 double omega = 0.0;
25 double phi = 0.0;
26
27 // 3D sphere
28
29 for (r = 0.3 ; r < 0.35 ;r += 0.001)
30 {
31 for (omega = 0.0; omega < M_PI ; omega += 0.006)
32 {
33 for (phi = 0.0; phi < 2.0*M_PI ; phi += 0.006)
34 {
35 Point<3,float> p;
36
37 p.get(0) = r*sin(omega)*sin(phi) + 0.5;
38 p.get(1) = r*sin(omega)*cos(phi) + 0.5;
39 p.get(2) = r*cos(omega) + 0.5;
40
41 // convert point into grid point
42
43 grid_key_dx<3> kd = cdsm.getCellGrid(p);
44
45 grid.template insert<0>(kd) = sin(omega)*sin(omega)*sin(2*phi);
46 grid.template insert<1>(kd) = 0;
47 }
48 }
49 }
50
51 auto it = grid.getIterator();
52
53 while (it.isNext())
54 {
55 tot_count++;
56
57 ++it;
58 }
59
60 return tot_count;
61}
62
63template <typename grid_type, typename cell_decomposer>
64size_t fill_sphere_quad(grid_type & grid, cell_decomposer & cdsm)
65{
66 size_t tot_count = 0;
67 double r = 0.3;
68 double omega = 0.0;
69 double phi = 0.0;
70
71 // 3D sphere
72
73 for (r = 0.3 ; r < 0.4 ;r += 0.001)
74 {
75 for (omega = 0.0; omega < M_PI ; omega += 0.006)
76 {
77 for (phi = 0.0; phi < 2.0*M_PI ; phi += 0.006)
78 {
79 Point<3,float> p;
80
81 p.get(0) = r*sin(omega)*sin(phi) + 0.5;
82 p.get(1) = r*sin(omega)*cos(phi) + 0.5;
83 p.get(2) = r*cos(omega) + 0.5;
84
85 // convert point into grid point
86
87 grid_key_dx<3> kd = cdsm.getCellGrid(p);
88
89 grid.template insert<0>(kd) = kd.get(0)*kd.get(0) + kd.get(1)*kd.get(1) + kd.get(2)*kd.get(2);
90 grid.template insert<1>(kd) = 0;
91 }
92 }
93 }
94
95 auto it = grid.getIterator();
96
97 while (it.isNext())
98 {
99 tot_count++;
100
101 ++it;
102 }
103
104 return tot_count;
105}
106
107template <typename grid_type, typename cell_decomposer>
108size_t fill_sphere_quad_v(grid_type & grid, cell_decomposer & cdsm)
109{
110 size_t tot_count = 0;
111 double r = 0.3;
112 double omega = 0.0;
113 double phi = 0.0;
114
115 // 3D sphere
116
117 for (r = 0.3 ; r < 0.4 ;r += 0.001)
118 {
119 for (omega = 0.0; omega < M_PI ; omega += 0.006)
120 {
121 for (phi = 0.0; phi < 2.0*M_PI ; phi += 0.006)
122 {
123 Point<3,float> p;
124
125 p.get(0) = r*sin(omega)*sin(phi) + 0.5;
126 p.get(1) = r*sin(omega)*cos(phi) + 0.5;
127 p.get(2) = r*cos(omega) + 0.5;
128
129 // convert point into grid point
130
131 grid_key_dx<3> kd = cdsm.getCellGrid(p);
132
133 grid.template insert<0>(kd) = kd.get(0)*kd.get(0) + kd.get(1)*kd.get(1) + kd.get(2)*kd.get(2);
134 grid.template insert<1>(kd) = 0;
135 grid.template insert<3>(kd)[0] = kd.get(0)*kd.get(0) + kd.get(1)*kd.get(1) + kd.get(2)*kd.get(2) + 10000;
136 grid.template insert<3>(kd)[1] = kd.get(0)*kd.get(0) + kd.get(1)*kd.get(1) + kd.get(2)*kd.get(2) + 60000;
137 grid.template insert<3>(kd)[2] = kd.get(0)*kd.get(0) + kd.get(1)*kd.get(1) + kd.get(2)*kd.get(2) + 80000;
138 }
139 }
140 }
141
142 auto it = grid.getIterator();
143
144 while (it.isNext())
145 {
146 tot_count++;
147
148 ++it;
149 }
150
151 return tot_count;
152}
153
154BOOST_AUTO_TEST_CASE( sparse_grid_use_test)
155{
156 size_t sz[3] = {10000,10000,10000};
157
158 sgrid_cpu<3,aggregate<float>,HeapMemory> grid(sz);
159
160 grid.getBackgroundValue().template get<0>() = 0.0;
161
162 // We fill a sphere with a band
163
164 grid_key_dx<3> key1({5000,5000,5000});
165 grid_key_dx<3> key2({5001,5001,5001});
166 grid_key_dx<3> key3({5002,5003,5003});
167
168 grid.template insert<0>(key1) = 1.0;
169 grid.template insert<0>(key2) = 2.0;
170 grid.template insert<0>(key3) = 3.0;
171
172 BOOST_REQUIRE_EQUAL(grid.template get<0>(key1),1.0);
173 BOOST_REQUIRE_EQUAL(grid.template get<0>(key2),2.0);
174 BOOST_REQUIRE_EQUAL(grid.template get<0>(key3),3.0);
175
176 auto it = grid.getIterator();
177
178 size_t count = 0;
179
180 while (it.isNext())
181 {
182 count++;
183
184 ++it;
185 }
186
187 BOOST_REQUIRE_EQUAL(count,(size_t)3);
188}
189
190BOOST_AUTO_TEST_CASE( sparse_grid_fill_all_test)
191{
192 size_t sz[3] = {171,171,171};
193
194 sgrid_cpu<3,aggregate<float>,HeapMemory> grid(sz);
195
196 grid.getBackgroundValue().template get<0>() = 0.0;
197
198 grid_sm<3,void> g_sm(sz);
199
200 grid_key_dx_iterator<3> kit(g_sm);
201
202 while (kit.isNext())
203 {
204 auto key = kit.get();
205
206 grid.template insert<0>(key) = g_sm.LinId(key);
207
208 ++kit;
209 }
210
211 auto it = grid.getIterator();
212
213 size_t count = 0;
214
215 bool match = true;
216
217 while (it.isNext())
218 {
219 auto key = it.get();
220
221 // return a grid_key_dx
222 auto key_pos = it.getKeyF();
223
224 match &= (grid.template get<0>(key_pos) == g_sm.LinId(key));
225
226 count++;
227
228 ++it;
229 }
230
231 BOOST_REQUIRE_EQUAL(count,(size_t)171*171*171);
232 BOOST_REQUIRE_EQUAL(grid.size(),(size_t)171*171*171);
233 BOOST_REQUIRE_EQUAL(match,true);
234
235 // remove all points
236
237 grid_key_dx_iterator<3> kit2(g_sm);
238
239 while (kit2.isNext())
240 {
241 auto key = kit2.get();
242
243 grid.remove(key);
244
245 ++kit2;
246 }
247
248 size_t tot = grid.size();
249 BOOST_REQUIRE_EQUAL(tot,0ul);
250}
251
252
253BOOST_AUTO_TEST_CASE( sparse_grid_fill_sparse_test)
254{
255 size_t sz[3] = {500,500,500};
256
257 sgrid_cpu<3,aggregate<double,int>,HeapMemory> grid(sz);
258
259 grid.getBackgroundValue().template get<0>() = 0.0;
260
261 CellDecomposer_sm<3, float, shift<3,float>> cdsm;
262
263 Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
264
265 cdsm.setDimensions(domain, sz, 0);
266
267 fill_sphere(grid,cdsm);
268
269 double r = 0.3;
270 double omega = 0.0;
271 double phi = 0.0;
272
273 for (r = 0.3 ; r < 0.35 ;r += 0.001)
274 {
275 for (omega = 0.0; omega < M_PI ; omega += 0.006)
276 {
277 for (phi = 0.0; phi < 2.0*M_PI ; phi += 0.006)
278 {
279
280 Point<3,float> p;
281
282 p.get(0) = r*sin(omega)*sin(phi) + 0.5;
283 p.get(1) = r*sin(omega)*cos(phi) + 0.5;
284 p.get(2) = r*cos(omega) + 0.5;
285
286 // convert point into grid point
287
288 grid_key_dx<3> kd = cdsm.getCellGrid(p);
289
290
291 if (grid.template get<0>(kd) == sin(omega)*sin(omega)*sin(2*phi))
292 {grid.template insert<1>(kd) = 1;}
293
294 }
295 }
296 }
297
298 auto it = grid.getIterator();
299
300 bool match = true;
301
302 while(it.isNext())
303 {
304 auto key = it.get();
305
306 if (grid.template get<1>(key) == 0)
307 {match = false;}
308
309 ++it;
310 }
311
312 BOOST_REQUIRE_EQUAL(match,true);
313
314 // remove the points
315
316 for (r = 0.3 ; r < 0.35 ;r += 0.001)
317 {
318 for (omega = 0.0; omega < M_PI ; omega += 0.006)
319 {
320 for (phi = 0.0; phi < 2.0*M_PI ; phi += 0.006)
321 {
322
323 Point<3,float> p;
324
325 p.get(0) = r*sin(omega)*sin(phi) + 0.5;
326 p.get(1) = r*sin(omega)*cos(phi) + 0.5;
327 p.get(2) = r*cos(omega) + 0.5;
328
329 // convert point into grid point
330
331 grid_key_dx<3> kd = cdsm.getCellGrid(p);
332
333
334 grid.remove(kd);
335
336 }
337 }
338 }
339
340 size_t tot;
341 tot = grid.size();
342
343 BOOST_REQUIRE_EQUAL(tot,0ul);
344}
345
346
347BOOST_AUTO_TEST_CASE( sparse_grid_resize_test)
348{
349 size_t sz[3] = {500,500,500};
350
351 sgrid_cpu<3,aggregate<double,int>,HeapMemory> grid(sz);
352 sgrid_cpu<3,aggregate<double,int>,HeapMemory> grid2(sz);
353
354 grid.getBackgroundValue().template get<0>() = 0.0;
355
356 CellDecomposer_sm<3, float, shift<3,float>> cdsm;
357
358 Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
359
360 cdsm.setDimensions(domain, sz, 0);
361
362 double r = 0.3;
363 double omega = 0.0;
364 double phi = 0.0;
365
366 // 3D sphere
367
368 for (r = 0.3 ; r < 0.35 ;r += 0.001)
369 {
370 for (omega = 0.0; omega < M_PI ; omega += 0.006)
371 {
372 for (phi = 0.0; phi < 2.0*M_PI ; phi += 0.006)
373 {
374 Point<3,float> p;
375
376 p.get(0) = r*sin(omega)*sin(phi) + 0.5;
377 p.get(1) = r*sin(omega)*cos(phi) + 0.5;
378 p.get(2) = r*cos(omega) + 0.5;
379
380 // convert point into grid point
381
382 grid_key_dx<3> kd = cdsm.getCellGrid(p);
383
384 grid.template insert<0>(kd) = sin(omega)*sin(omega)*sin(2*phi);
385 grid.template insert<1>(kd) = 0;
386 grid2.template insert<0>(kd) = sin(omega)*sin(omega)*sin(2*phi);
387 grid2.template insert<1>(kd) = 0;
388 }
389 }
390 }
391
392 size_t sz_b[3] = {1024,1024,1024};
393 grid2.resize(sz_b);
394
395 // Check that both grid contain the same information
396
397 auto it = grid2.getIterator();
398
399 bool match = true;
400
401 while(it.isNext())
402 {
403 auto key = it.get();
404
405 if (grid.template get<0>(key) != grid2.template get<0>(key))
406 {
407 match = false;
408 break;
409 }
410
411 ++it;
412 }
413
414 BOOST_REQUIRE_EQUAL(match,true);
415
416 // now we resize smalle
417
418 size_t sz_s[3] = {250,250,250};
419 grid2.resize(sz_s);
420
421 //
422
423 auto it2 = grid.getIterator();
424
425 match = true;
426
427 while(it2.isNext())
428 {
429 auto key = it2.get();
430
431 // we check if the key is inside
432
433 bool cin = true;
434
435 cin &= (size_t)key.get(0) < sz_s[0];
436 cin &= (size_t)key.get(1) < sz_s[1];
437 cin &= (size_t)key.get(2) < sz_s[2];
438
439
440 if (cin == true)
441 {
442 if (grid.template get<0>(key) != grid2.template get<0>(key))
443 {
444 match = false;
445 break;
446 }
447 }
448
449 ++it2;
450 }
451
452 BOOST_REQUIRE_EQUAL(match,true);
453}
454
455
456
457BOOST_AUTO_TEST_CASE( sparse_grid_fill_all_with_resize_test)
458{
459 size_t sz[3] = {10,10,171};
460
461 sgrid_cpu<3,aggregate<float>,HeapMemory> grid(sz);
462
463 grid.getBackgroundValue().template get<0>() = 0.0;
464
465 grid_sm<3,void> g_sm(sz);
466
467 grid_key_dx_iterator<3> kit(g_sm);
468
469 while (kit.isNext())
470 {
471 auto key = kit.get();
472
473 grid.template insert<0>(key) = g_sm.LinId(key);
474
475 ++kit;
476 }
477
478 size_t sz_b[3] = {20,20,200};
479
480 // now we increase the size
481 grid.resize(sz_b);
482
483 auto it = grid.getIterator();
484
485 size_t count = 0;
486
487 bool match = true;
488
489 while (it.isNext())
490 {
491 auto key = it.get();
492
493 // return a grid_key_dx
494 auto key_pos = it.getKeyF();
495
496 match &= (grid.template get<0>(key_pos) == g_sm.LinId(key));
497
498 count++;
499
500 ++it;
501 }
502
503 BOOST_REQUIRE_EQUAL(count,(size_t)10*10*171);
504 BOOST_REQUIRE_EQUAL(grid.size(),(size_t)10*10*171);
505 BOOST_REQUIRE_EQUAL(match,true);
506
507 // refill with the full set of point
508
509 grid_sm<3,void> g_sm2(sz_b);
510
511 grid_key_dx_iterator<3> kit2(g_sm2);
512
513 while (kit2.isNext())
514 {
515 auto key = kit2.get();
516
517 grid.template insert<0>(key) = g_sm2.LinId(key);
518
519 ++kit2;
520 }
521
522 auto it2 = grid.getIterator();
523
524 count = 0;
525
526 match = true;
527
528 while (it2.isNext())
529 {
530 auto key = it2.get();
531
532 // return a grid_key_dx
533 auto key_pos = it2.getKeyF();
534
535 match &= (grid.template get<0>(key_pos) == g_sm2.LinId(key));
536
537 count++;
538
539 ++it2;
540 }
541
542 BOOST_REQUIRE_EQUAL(count,(size_t)20*20*200);
543 BOOST_REQUIRE_EQUAL(grid.size(),(size_t)20*20*200);
544 BOOST_REQUIRE_EQUAL(match,true);
545}
546
547BOOST_AUTO_TEST_CASE( sparse_grid_insert_o_test)
548{
549 size_t sz[3] = {10,10,171};
550
551 sgrid_cpu<3,aggregate<float>,HeapMemory> grid(sz);
552
553 grid.getBackgroundValue().template get<0>() = 0.0;
554
555 size_t ele;
556 grid_key_dx<3> key({5,5,90});
557
558 auto & flt = grid.insert_o(key,ele).template get<0>();
559 flt[ele] = 117.0;
560
561 BOOST_REQUIRE_EQUAL(grid.template get<0>(key),117.0);
562}
563
564
565BOOST_AUTO_TEST_CASE( sparse_grid_sub_grid_it)
566{
567 size_t sz[3] = {171,171,171};
568
569 sgrid_cpu<3,aggregate<float>,HeapMemory> grid(sz);
570
571 grid.getBackgroundValue().template get<0>() = 0.0;
572
573 grid_sm<3,void> g_sm(sz);
574
575 grid_key_dx_iterator<3> kit(g_sm);
576
577 while (kit.isNext())
578 {
579 auto key = kit.get();
580
581 grid.template insert<0>(key) = g_sm.LinId(key);
582
583 ++kit;
584 }
585
586 grid_key_dx<3> start({21,21,21});
587 grid_key_dx<3> stop({90,90,90});
588
589 bool error = false;
590 size_t count = 0;
591 auto it_sub = grid.getIterator(start,stop);
592
593 while (it_sub.isNext())
594 {
595 auto gkey = it_sub.get();
596
597 if (gkey.get(0) < start.get(0) ||
598 gkey.get(1) < start.get(1) ||
599 gkey.get(2) < start.get(2) ||
600 gkey.get(0) > stop.get(0) ||
601 gkey.get(1) > stop.get(1) ||
602 gkey.get(2) > stop.get(2))
603 {
604 error = true;
605 }
606
607 count++;
608
609 ++it_sub;
610 }
611
612 size_t tot = (stop.get(2) - start.get(2) + 1)*(stop.get(1) - start.get(1) + 1)*(stop.get(0) - start.get(0) + 1);
613 BOOST_REQUIRE_EQUAL(error,false);
614 BOOST_REQUIRE_EQUAL(count,tot);
615}
616
617
618BOOST_AUTO_TEST_CASE( sparse_grid_sub_grid_it_quarter_sphere)
619{
620 size_t sz[3] = {501,501,501};
621 size_t sz_cell[3] = {500,500,500};
622
623 sgrid_cpu<3,aggregate<double,int>,HeapMemory> grid(sz);
624
625 grid.getBackgroundValue().template get<0>() = 0.0;
626
627 CellDecomposer_sm<3, float, shift<3,float>> cdsm;
628
629 Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
630
631 cdsm.setDimensions(domain, sz_cell, 0);
632
633 fill_sphere(grid,cdsm);
634
635 grid_key_dx<3> start({0,0,0});
636 grid_key_dx<3> stop({250,250,250});
637
638 bool error = false;
639 size_t count = 0;
640 auto it_sub = grid.getIterator(start,stop);
641
642 while (it_sub.isNext())
643 {
644 auto gkey = it_sub.get();
645
646 if (gkey.get(0) < start.get(0) ||
647 gkey.get(1) < start.get(1) ||
648 gkey.get(2) < start.get(2) ||
649 gkey.get(0) > stop.get(0) ||
650 gkey.get(1) > stop.get(1) ||
651 gkey.get(2) > stop.get(2))
652 {
653 error = true;
654 }
655
656 // Check that the point is in the sphere
657
658 double radius = (gkey.get(0) - 250)*(gkey.get(0) - 250) +
659 (gkey.get(1) - 250)*(gkey.get(1) - 250) +
660 (gkey.get(2) - 250)*(gkey.get(2) - 250);
661
662 radius = sqrt(radius);
663
664 if (radius < 150 || radius >= 175)
665 {
666 // if is not in the radius remove it
667 grid.remove(gkey);
668 }
669
670 count++;
671
672 ++it_sub;
673 }
674
675 BOOST_REQUIRE_EQUAL(error,false);
676
677 // We go again across the point now every point out the sphere is an error
678
679 count = 0;
680 auto it_sub2 = grid.getIterator(start,stop);
681
682 while (it_sub2.isNext())
683 {
684 auto gkey = it_sub2.get();
685
686 if (gkey.get(0) < start.get(0) ||
687 gkey.get(1) < start.get(1) ||
688 gkey.get(2) < start.get(2) ||
689 gkey.get(0) > stop.get(0) ||
690 gkey.get(1) > stop.get(1) ||
691 gkey.get(2) > stop.get(2))
692 {
693 error = true;
694 }
695
696 // Check that the point is in the sphere
697
698 double radius = (gkey.get(0) - 250)*(gkey.get(0) - 250) +
699 (gkey.get(1) - 250)*(gkey.get(1) - 250) +
700 (gkey.get(2) - 250)*(gkey.get(2) - 250);
701
702 radius = sqrt(radius);
703
704 if (radius < 150 || radius >= 175)
705 {
706 error = true;
707 }
708
709 count++;
710
711 ++it_sub2;
712 }
713
714 BOOST_REQUIRE_EQUAL(error,false);
715}
716
717BOOST_AUTO_TEST_CASE( sparse_grid_fast_stencil)
718{
719 size_t sz[3] = {501,501,501};
720 size_t sz_cell[3] = {500,500,500};
721
722 sgrid_soa<3,aggregate<double,double,int>,HeapMemory> grid(sz);
723
724 grid.getBackgroundValue().template get<0>() = 0.0;
725
726 CellDecomposer_sm<3, float, shift<3,float>> cdsm;
727
728 Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
729
730 cdsm.setDimensions(domain, sz_cell, 0);
731
732 fill_sphere_quad(grid,cdsm);
733
734 grid_key_dx<3> start({1,1,1});
735 grid_key_dx<3> stop({499,499,499});
736
737 for (int i = 0 ; i < 1 ; i++)
738 {
739
740 grid.private_get_nnlist().resize(NNStar_c<3>::nNN * grid.private_get_header_inf().size());
741 auto it = grid.getBlockIterator<1>(start,stop);
742
743 unsigned char mask[decltype(it)::sizeBlockBord];
744 __attribute__ ((aligned (32))) double block_bord_src[decltype(it)::sizeBlockBord];
745 __attribute__ ((aligned (32))) double block_bord_dst[decltype(it)::sizeBlock];
746
747 const Vc::double_v six(6.0);
748
749 while (it.isNext())
750 {
751 it.loadBlockBorder<0,NNStar_c<3>,false>(block_bord_src,mask);
752
753 for (int i = 0 ; i < 1 ; i++)
754 {
755
756 for (int k = it.start_b(2) ; k < it.stop_b(2) ; k++)
757 {
758 for (int j = it.start_b(1) ; j < it.stop_b(1) ; j++)
759 {
760 int c = it.LinB(it.start_b(0),j,k);
761
762 int xp = c+1;
763 int xm = c-1;
764
765 int yp = it.LinB(it.start_b(0),j+1,k);
766 int ym = it.LinB(it.start_b(0),j-1,k);
767
768 int zp = it.LinB(it.start_b(0),j,k+1);
769 int zm = it.LinB(it.start_b(0),j,k-1);
770
771 for (int i = it.start_b(0) ; i < it.stop_b(0) ; i++)
772 {
773 // we do only id exist the point
774 if (mask[c] == false) {continue;}
775
776 bool surround = mask[xp] & mask[xm] & mask[ym] & mask[yp] & mask[zp] & mask[zm];
777
778 double Lap = block_bord_src[xp] + block_bord_src[xm] +
779 block_bord_src[yp] + block_bord_src[ym] +
780 block_bord_src[zp] + block_bord_src[zm] - 6.0*block_bord_src[c];
781
782 Lap = (surround)?Lap:6.0;
783
784 block_bord_dst[it.LinB_off(i,j,k)] = Lap;
785
786 c++;
787 xp++;
788 xm++;
789 yp++;
790 ym++;
791 zp++;
792 zm++;
793 }
794 }
795 }
796
797 }
798
799 it.storeBlock<1>(block_bord_dst);
800
801 ++it;
802 }
803
804 }
805
806 bool check = true;
807 auto it2 = grid.getIterator();
808 while (it2.isNext())
809 {
810 auto p = it2.get();
811
812 check &= grid.template get<1>(p) == 6;
813
814 ++it2;
815 }
816
817 BOOST_REQUIRE_EQUAL(check,true);
818 // Check correct-ness
819
820// print_grid("debug_out",grid);
821}
822
823BOOST_AUTO_TEST_CASE( sparse_grid_fast_stencil_vectorized)
824{
825 size_t sz[3] = {501,501,501};
826 size_t sz_cell[3] = {500,500,500};
827
828 sgrid_soa<3,aggregate<double,double,int>,HeapMemory> grid(sz);
829
830 grid.getBackgroundValue().template get<0>() = 0.0;
831
832 CellDecomposer_sm<3, float, shift<3,float>> cdsm;
833
834 Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
835
836 cdsm.setDimensions(domain, sz_cell, 0);
837
838 fill_sphere_quad(grid,cdsm);
839
840 grid_key_dx<3> start({1,1,1});
841 grid_key_dx<3> stop({499,499,499});
842
843 grid.private_get_nnlist().resize(NNStar_c<3>::nNN * grid.private_get_header_inf().size());
844
845 for (int i = 0 ; i < 1 ; i++)
846 {
847
848 auto it = grid.getBlockIterator<1>(start,stop);
849
850 unsigned char mask[decltype(it)::sizeBlockBord];
851 unsigned char mask_sum[decltype(it)::sizeBlockBord];
852 __attribute__ ((aligned (32))) double block_bord_src[decltype(it)::sizeBlockBord];
853 __attribute__ ((aligned (32))) double block_bord_dst[decltype(it)::sizeBlock];
854
855 typedef typename boost::mpl::at<decltype(it)::stop_border_vmpl,boost::mpl::int_<0>>::type sz0;
856 typedef typename boost::mpl::at<decltype(it)::stop_border_vmpl,boost::mpl::int_<1>>::type sz1;
857 typedef typename boost::mpl::at<decltype(it)::stop_border_vmpl,boost::mpl::int_<2>>::type sz2;
858
859 const Vc::double_v six(6.0);
860
861 while (it.isNext())
862 {
863 it.loadBlockBorder<0,NNStar_c<3>,false>(block_bord_src,mask);
864
865 // Sum the mask
866 for (int k = it.start_b(2) ; k < it.stop_b(2) ; k++)
867 {
868 for (int j = it.start_b(1) ; j < it.stop_b(1) ; j++)
869 {
870 int c = it.LinB(it.start_b(0),j,k);
871
872 int yp = it.LinB(it.start_b(0),j+1,k);
873 int ym = it.LinB(it.start_b(0),j-1,k);
874
875 int zp = it.LinB(it.start_b(0),j,k+1);
876 int zm = it.LinB(it.start_b(0),j,k-1);
877
878 for (int i = it.start_b(0) ; i < it.stop_b(0) ; i += sizeof(size_t))
879 {
880 size_t cmd = *(size_t *)&mask[c];
881
882 if (cmd == 0) {continue;}
883
884 size_t xpd = *(size_t *)&mask[c+1];
885 size_t xmd = *(size_t *)&mask[c-1];
886 size_t ypd = *(size_t *)&mask[yp];
887 size_t ymd = *(size_t *)&mask[ym];
888 size_t zpd = *(size_t *)&mask[zp];
889 size_t zmd = *(size_t *)&mask[zm];
890
891 size_t sum = xpd + xmd +
892 ypd + ymd +
893 zpd + zmd + cmd;
894
895 *(size_t *)&mask_sum[c] = sum;
896
897 c += sizeof(size_t);
898 yp += sizeof(size_t);
899 ym += sizeof(size_t);
900 zp += sizeof(size_t);
901 zm += sizeof(size_t);
902 }
903 }
904 }
905
906 for (int k = it.start_b(2) ; k < it.stop_b(2) ; k++)
907 {
908 for (int j = it.start_b(1) ; j < it.stop_b(1) ; j++)
909 {
910 int c = it.LinB(it.start_b(0),j,k);
911
912 int yp = it.LinB(it.start_b(0),j+1,k);
913 int ym = it.LinB(it.start_b(0),j-1,k);
914
915 int zp = it.LinB(it.start_b(0),j,k+1);
916 int zm = it.LinB(it.start_b(0),j,k-1);
917
918 int cd = it.LinB_off(it.start_b(0),j,k);
919
920 for (int i = it.start_b(0) ; i < it.stop_b(0) ; i += Vc::double_v::Size)
921 {
922 Vc::Mask<double> cmp;
923
924 if (Vc::double_v::Size == 2)
925 {
926 cmp[0] = mask[c] == true;
927 cmp[1] = mask[c+1] == true;
928 }
929 else if (Vc::double_v::Size == 4)
930 {
931 cmp[0] = mask[c] == true;
932 cmp[1] = mask[c+1] == true;
933 cmp[2] = mask[c+2] == true;
934 cmp[3] = mask[c+3] == true;
935 }
936 else
937 {
938 std::cout << "UNSUPPORTED" << std::endl;
939 exit(1);
940 }
941
942
943 // we do only id exist the point
944 if (Vc::none_of(cmp) == true) {continue;}
945
946 Vc::Mask<double> surround;
947
948 Vc::double_v xpd(&block_bord_src[c+1],Vc::Unaligned);
949 Vc::double_v xmd(&block_bord_src[c-1],Vc::Unaligned);
950 Vc::double_v ypd(&block_bord_src[c+sz0::value],Vc::Unaligned);
951 Vc::double_v ymd(&block_bord_src[c-sz0::value],Vc::Unaligned);
952 Vc::double_v zpd(&block_bord_src[zp],Vc::Unaligned);
953 Vc::double_v zmd(&block_bord_src[zm],Vc::Unaligned);
954 Vc::double_v cmd(&block_bord_src[c],Vc::Unaligned);
955
956 Vc::double_v Lap = xpd + xmd +
957 ypd + ymd +
958 zpd + zmd - 6.0*cmd;
959
960 if (Vc::double_v::Size == 2)
961 {
962 surround[0] = (mask_sum[c] == 7);
963 surround[1] = (mask_sum[c+1] == 7);
964 }
965 else if (Vc::double_v::Size == 4)
966 {
967 surround[0] = (mask_sum[c] == 7);
968 surround[1] = (mask_sum[c+1] == 7);
969 surround[2] = (mask_sum[c+2] == 7);
970 surround[3] = (mask_sum[c+3] == 7);
971 }
972
973 Lap = Vc::iif(surround,Lap,six);
974
975 Lap.store(&block_bord_dst[cd],cmp,Vc::Aligned);
976
977 c+=Vc::double_v::Size;
978 zp+=Vc::double_v::Size;
979 zm+=Vc::double_v::Size;
980 cd+=Vc::double_v::Size;
981 }
982 }
983 }
984
985 it.storeBlock<1>(block_bord_dst);
986
987 ++it;
988 }
989
990 }
991
992 bool check = true;
993 auto it2 = grid.getIterator();
994 while (it2.isNext())
995 {
996 auto p = it2.get();
997
998 check &= grid.template get<1>(p) == 6;
999
1000 ++it2;
1001 }
1002
1003 BOOST_REQUIRE_EQUAL(check,true);
1004 // Check correct-ness
1005
1006 //print_grid("debug_out",grid);
1007}
1008
1009BOOST_AUTO_TEST_CASE( sparse_grid_fast_stencil_vectorized_block_skip)
1010{
1011 size_t sz[3] = {501,501,501};
1012 size_t sz_cell[3] = {500,500,500};
1013
1014 sgrid_soa<3,aggregate<double,double,int>,HeapMemory> grid(sz);
1015
1016 grid.getBackgroundValue().template get<0>() = 0.0;
1017
1018 CellDecomposer_sm<3, float, shift<3,float>> cdsm;
1019
1020 Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
1021
1022 cdsm.setDimensions(domain, sz_cell, 0);
1023
1024 fill_sphere_quad(grid,cdsm);
1025
1026 grid.reorder();
1027
1028 grid_key_dx<3> start({1,1,1});
1029 grid_key_dx<3> stop({499,499,499});
1030
1031 for (int i = 0 ; i < 1 ; i++)
1032 {
1033
1034 auto it = grid.getBlockIterator<1>(start,stop);
1035
1036 auto & datas = grid.private_get_data();
1037 auto & headers = grid.private_get_header_mask();
1038
1039 typedef typename boost::mpl::at<decltype(it)::stop_border_vmpl,boost::mpl::int_<0>>::type sz0;
1040 typedef typename boost::mpl::at<decltype(it)::stop_border_vmpl,boost::mpl::int_<1>>::type sz1;
1041 typedef typename boost::mpl::at<decltype(it)::stop_border_vmpl,boost::mpl::int_<2>>::type sz2;
1042
1043 typedef typename sgrid_soa<3,aggregate<double,double,int>,HeapMemory>::chunking_type chunking;
1044
1045 const Vc::double_v one(1.0);
1046
1047 while (it.isNext())
1048 {
1049 // Load
1050 long int offset_jump[6];
1051 size_t cid = it.getChunkId();
1052
1053 auto chunk = datas.get(cid);
1054 auto & mask = headers.get(cid);
1055
1056 bool exist;
1057 grid_key_dx<3> p = grid.getChunkPos(cid) + grid_key_dx<3>({-1,0,0});
1058 long int r = grid.getChunk(p,exist);
1059 offset_jump[0] = (r-cid)*decltype(it)::sizeBlock;
1060
1061 p = grid.getChunkPos(cid) + grid_key_dx<3>({1,0,0});
1062 r = grid.getChunk(p,exist);
1063 offset_jump[1] = (r-cid)*decltype(it)::sizeBlock;
1064
1065 p = grid.getChunkPos(cid) + grid_key_dx<3>({0,-1,0});
1066 r = grid.getChunk(p,exist);
1067 offset_jump[2] = (r-cid)*decltype(it)::sizeBlock;
1068
1069 p = grid.getChunkPos(cid) + grid_key_dx<3>({0,1,0});
1070 r = grid.getChunk(p,exist);
1071 offset_jump[3] = (r-cid)*decltype(it)::sizeBlock;
1072
1073 p = grid.getChunkPos(cid) + grid_key_dx<3>({0,0,-1});
1074 r = grid.getChunk(p,exist);
1075 offset_jump[4] = (r-cid)*decltype(it)::sizeBlock;
1076
1077 p = grid.getChunkPos(cid) + grid_key_dx<3>({0,0,1});
1078 r = grid.getChunk(p,exist);
1079 offset_jump[5] = (r-cid)*decltype(it)::sizeBlock;
1080
1081 // Load offset jumps
1082
1083 long int s2 = 0;
1084
1085 typedef boost::mpl::at<typename chunking::type,boost::mpl::int_<2>>::type sz;
1086 typedef boost::mpl::at<typename chunking::type,boost::mpl::int_<1>>::type sy;
1087 typedef boost::mpl::at<typename chunking::type,boost::mpl::int_<0>>::type sx;
1088
1089 for (int v = 0 ; v < sz::value ; v++)
1090 {
1091 for (int j = 0 ; j < sy::value ; j++)
1092 {
1093 for (int k = 0 ; k < sx::value ; k += Vc::double_v::Size)
1094 {
1095 // we do only id exist the point
1096 if (*(int *)&mask.mask[s2] == 0) {s2 += Vc::double_v::Size; continue;}
1097
1098 Vc::Mask<double> surround;
1099
1100 data_il<Vc::double_v::Size> mxm;
1101 data_il<Vc::double_v::Size> mxp;
1102 data_il<Vc::double_v::Size> mym;
1103 data_il<Vc::double_v::Size> myp;
1104 data_il<Vc::double_v::Size> mzm;
1105 data_il<Vc::double_v::Size> mzp;
1106
1107 Vc::double_v xm;
1108 Vc::double_v xp;
1109
1110 Vc::double_v cmd(&chunk.template get<0>()[s2]);
1111
1112 // Load x-1
1113 long int sumxm = s2-1;
1114 sumxm += (k==0)?offset_jump[0] + sx::value:0;
1115
1116 // Load x+1
1117 long int sumxp = s2+Vc::double_v::Size;
1118 sumxp += (k+Vc::double_v::Size == sx::value)?offset_jump[1] - sx::value:0;
1119
1120 long int sumym = (j == 0)?offset_jump[2] + (sy::value-1)*sx::value:-sx::value;
1121 sumym += s2;
1122 long int sumyp = (j == sy::value-1)?offset_jump[3] - (sy::value - 1)*sx::value:sx::value;
1123 sumyp += s2;
1124 long int sumzm = (v == 0)?offset_jump[4] + (sz::value-1)*sx::value*sy::value:-sx::value*sy::value;
1125 sumzm += s2;
1126 long int sumzp = (v == sz::value-1)?offset_jump[5] - (sz::value - 1)*sx::value*sy::value:sx::value*sy::value;
1127 sumzp += s2;
1128
1129 if (Vc::double_v::Size == 2)
1130 {
1131 mxm.uc[0] = mask.mask[sumxm];
1132 mxm.uc[1] = mask.mask[s2];
1133
1134 mxp.uc[0] = mask.mask[s2+1];
1135 mxp.uc[1] = mask.mask[sumxp];
1136 }
1137 else if (Vc::double_v::Size == 4)
1138 {
1139 mxm.uc[0] = mask.mask[sumxm];
1140 mxm.uc[1] = mask.mask[s2];
1141 mxm.uc[2] = mask.mask[s2+1];
1142 mxm.uc[3] = mask.mask[s2+2];
1143
1144 mxp.uc[0] = mask.mask[s2+1];
1145 mxp.uc[1] = mask.mask[s2+2];
1146 mxp.uc[2] = mask.mask[s2+3];
1147 mxp.uc[3] = mask.mask[sumxp];
1148 }
1149 else
1150 {
1151 std::cout << "UNSUPPORTED" << std::endl;
1152 exit(1);
1153 }
1154
1155 mym.i = *(int *)&mask.mask[sumym];
1156 myp.i = *(int *)&mask.mask[sumyp];
1157
1158 mzm.i = *(int *)&mask.mask[sumzm];
1159 mzp.i = *(int *)&mask.mask[sumzp];
1160
1161 if (Vc::double_v::Size == 2)
1162 {
1163 xm[0] = chunk.template get<0>()[sumxm];
1164 xm[1] = cmd[0];
1165
1166 xp[0] = cmd[1];
1167 xp[1] = chunk.template get<0>()[sumxp];
1168 }
1169 else if (Vc::double_v::Size == 4)
1170 {
1171 xm[0] = chunk.template get<0>()[sumxm];
1172 xm[1] = cmd[0];
1173 xm[2] = cmd[1];
1174 xm[3] = cmd[2];
1175
1176 xp[0] = cmd[1];
1177 xp[1] = cmd[2];
1178 xp[2] = cmd[3];
1179 xp[3] = chunk.template get<0>()[sumxp];
1180 }
1181
1182 // Load y and z direction
1183
1184 Vc::double_v ym(&chunk.template get<0>()[sumym],Vc::Aligned);
1185 Vc::double_v yp(&chunk.template get<0>()[sumyp],Vc::Aligned);
1186 Vc::double_v zm(&chunk.template get<0>()[sumzm],Vc::Aligned);
1187 Vc::double_v zp(&chunk.template get<0>()[sumzp],Vc::Aligned);
1188
1189 // Calculate
1190
1191 data_il<Vc::double_v::Size> tot_m;
1192 tot_m.i = mxm.i + mxp.i + mym.i + myp.i + mzm.i + mzp.i;
1193
1194 if (Vc::double_v::Size == 2)
1195 {
1196 surround[0] = (tot_m.uc[0] == 6);
1197 surround[1] = (tot_m.uc[1] == 6);
1198 }
1199 else if (Vc::double_v::Size == 4)
1200 {
1201 surround[0] = (tot_m.uc[0] == 6);
1202 surround[1] = (tot_m.uc[1] == 6);
1203 surround[2] = (tot_m.uc[2] == 6);
1204 surround[3] = (tot_m.uc[3] == 6);
1205 }
1206
1207 Vc::double_v Lap = xp + xm +
1208 yp + ym +
1209 zp + zm - 6.0*cmd;
1210
1211 Lap = Vc::iif(surround,Lap,one);
1212
1213 Lap.store(&chunk.template get<1>()[s2],Vc::Aligned);
1214
1215 s2 += Vc::double_v::Size;
1216 }
1217 }
1218 }
1219
1220 ++it;
1221 }
1222
1223 }
1224
1225 int tot_six = 0;
1226 int tot_one = 0;
1227
1228 bool check = true;
1229 auto it2 = grid.getIterator();
1230 while (it2.isNext())
1231 {
1232 auto p = it2.get();
1233
1234 check &= (grid.template get<1>(p) == 6 || grid.template get<1>(p) == 1);
1235
1236 if (grid.template get<1>(p) == 1)
1237 {
1238 tot_one++;
1239 }
1240
1241 if (grid.template get<1>(p) == 6)
1242 {
1243 tot_six++;
1244 }
1245
1246 ++it2;
1247 }
1248
1249 BOOST_REQUIRE_EQUAL(check,true);
1250 BOOST_REQUIRE_EQUAL(tot_six,15857813);
1251 BOOST_REQUIRE_EQUAL(tot_one,2977262);
1252 // Check correct-ness
1253
1254 //print_grid("debug_out",grid);
1255}
1256
1257BOOST_AUTO_TEST_CASE( sparse_grid_fast_stencil_vectorized_simplified)
1258{
1259 size_t sz[3] = {501,501,501};
1260 size_t sz_cell[3] = {500,500,500};
1261
1262 sgrid_soa<3,aggregate<double,double,int>,HeapMemory> grid(sz);
1263
1264 grid.getBackgroundValue().template get<0>() = 0.0;
1265
1266 CellDecomposer_sm<3, float, shift<3,float>> cdsm;
1267
1268 Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
1269
1270 cdsm.setDimensions(domain, sz_cell, 0);
1271
1272 fill_sphere_quad(grid,cdsm);
1273
1274 grid_key_dx<3> start({1,1,1});
1275 grid_key_dx<3> stop({250,250,250});
1276
1277 for (int i = 0 ; i < 1 ; i++)
1278 {
1279
1280 int stencil[6][3] = {{1,0,0},{-1,0,0},{0,-1,0},{0,1,0},{0,0,-1},{0,0,1}};
1281
1282 grid.conv<0,1,1>(stencil,start,stop,[](Vc::double_v (& xs)[7], unsigned char * mask_sum){
1283 Vc::double_v Lap = xs[1] + xs[2] +
1284 xs[3] + xs[4] +
1285 xs[5] + xs[6] - 6.0*xs[0];
1286
1287 auto surround = load_mask<Vc::double_v>(mask_sum);
1288
1289 auto sur_mask = (surround == 6.0);
1290
1291 Lap = Vc::iif(sur_mask,Lap,Vc::double_v(1.0));
1292
1293 return Lap;
1294 });
1295 }
1296
1297 int tot_six = 0;
1298 int tot_one = 0;
1299
1300 bool check = true;
1301 auto it2 = grid.getIterator(start,stop);
1302 while (it2.isNext())
1303 {
1304 auto p = it2.get();
1305
1306 check &= (grid.template get<1>(p) == 6 || grid.template get<1>(p) == 1);
1307
1308 // Check the six should be a six
1309 auto xp = p.move(0,1);
1310 auto xm = p.move(0,-1);
1311
1312 auto yp = p.move(1,1);
1313 auto ym = p.move(1,-1);
1314
1315 auto zp = p.move(2,1);
1316 auto zm = p.move(2,-1);
1317
1318 bool is_six;
1319 if (grid.existPoint(xp) && grid.existPoint(xm) &&
1320 grid.existPoint(yp) && grid.existPoint(ym) &&
1321 grid.existPoint(zp) && grid.existPoint(zm))
1322 {
1323 is_six = true;
1324 }
1325 else
1326 {
1327 is_six = false;
1328 }
1329
1330 if (is_six == true && grid.template get<1>(p) != 6.0)
1331 {
1332 check = false;
1333 break;
1334 }
1335
1336 if (grid.template get<1>(p) == 1)
1337 {
1338 tot_one++;
1339 }
1340
1341 if (grid.template get<1>(p) == 6)
1342 {
1343 tot_six++;
1344 }
1345
1346 ++it2;
1347 }
1348
1349 BOOST_REQUIRE_EQUAL(check,true);
1350 BOOST_REQUIRE_EQUAL(tot_six,2020091);
1351 BOOST_REQUIRE_EQUAL(tot_one,376791);
1352 // Check correct-ness
1353
1354// print_grid("debug_out.vtk",grid);
1355}
1356
1357BOOST_AUTO_TEST_CASE( sparse_grid_fast_stencil_vectorized_cross_simplified)
1358{
1359 size_t sz[3] = {501,501,501};
1360 size_t sz_cell[3] = {500,500,500};
1361
1362 sgrid_soa<3,aggregate<double,double,int>,HeapMemory> grid(sz);
1363
1364 grid.getBackgroundValue().template get<0>() = 0.0;
1365
1366 CellDecomposer_sm<3, float, shift<3,float>> cdsm;
1367
1368 Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
1369
1370 cdsm.setDimensions(domain, sz_cell, 0);
1371
1372 fill_sphere_quad(grid,cdsm);
1373
1374 //grid.reorder();
1375
1376 grid_key_dx<3> start({1,1,1});
1377 grid_key_dx<3> stop({499,499,499});
1378
1379 for (int i = 0 ; i < 1 ; i++)
1380 {
1381 grid.conv_cross<0,1,1>(start,stop,[]( Vc::double_v & cmd, cross_stencil_v & s,
1382 unsigned char * mask_sum){
1383
1384 Vc::double_v Lap = s.xm + s.xp +
1385 s.ym + s.yp +
1386 s.zm + s.zp - 6.0*cmd;
1387
1388 Vc::Mask<double> surround;
1389
1390 for (int i = 0 ; i < Vc::double_v::Size ; i++)
1391 {surround[i] = (mask_sum[i] == 6);}
1392
1393 Lap = Vc::iif(surround,Lap,Vc::double_v(1.0));
1394
1395 return Lap;
1396 });
1397 }
1398
1399 int tot_six = 0;
1400 int tot_one = 0;
1401
1402 bool check = true;
1403 auto it2 = grid.getIterator(start,stop);
1404 while (it2.isNext())
1405 {
1406 auto p = it2.get();
1407
1408 check &= (grid.template get<1>(p) == 6 || grid.template get<1>(p) == 1);
1409
1410 // Check the six should be a six
1411 auto xp = p.move(0,1);
1412 auto xm = p.move(0,-1);
1413
1414 auto yp = p.move(1,1);
1415 auto ym = p.move(1,-1);
1416
1417 auto zp = p.move(2,1);
1418 auto zm = p.move(2,-1);
1419
1420 bool is_six;
1421 if (grid.existPoint(xp) && grid.existPoint(xm) &&
1422 grid.existPoint(yp) && grid.existPoint(ym) &&
1423 grid.existPoint(zp) && grid.existPoint(zm))
1424 {
1425 is_six = true;
1426 }
1427 else
1428 {
1429 is_six = false;
1430 }
1431
1432 if (is_six == true && grid.template get<1>(p) != 6.0)
1433 {
1434 check = false;
1435 break;
1436 }
1437
1438 if (grid.template get<1>(p) == 1)
1439 {
1440 tot_one++;
1441 }
1442
1443 if (grid.template get<1>(p) == 6)
1444 {
1445 tot_six++;
1446 }
1447
1448 ++it2;
1449 }
1450
1451 BOOST_REQUIRE_EQUAL(check,true);
1452 BOOST_REQUIRE_EQUAL(tot_six,15857813);
1453 BOOST_REQUIRE_EQUAL(tot_one,2977262);
1454 // Check correct-ness
1455
1456// print_grid("debug_out",grid);
1457}
1458
1459constexpr int x = 0;
1460constexpr int y = 1;
1461constexpr int z = 2;
1462
1463BOOST_AUTO_TEST_CASE( sparse_grid_fast_stencil_vectorized_cross_simplified_ids)
1464{
1465 size_t sz[3] = {501,501,501};
1466 size_t sz_cell[3] = {500,500,500};
1467
1468 sgrid_soa<3,aggregate<double,double,int>,HeapMemory> grid(sz);
1469
1470 grid.getBackgroundValue().template get<0>() = 0.0;
1471
1472 CellDecomposer_sm<3, float, shift<3,float>> cdsm;
1473
1474 Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
1475
1476 cdsm.setDimensions(domain, sz_cell, 0);
1477
1478 fill_sphere_quad(grid,cdsm);
1479
1480 //grid.reorder();
1481
1482 grid_key_dx<3> start({1,1,1});
1483 grid_key_dx<3> stop({499,499,499});
1484
1485 for (int i = 0 ; i < 1 ; i++)
1486 {
1487 grid.conv_cross_ids<1,double>(start,stop,[](auto & grid, auto & ids,
1488 unsigned char * mask_sum){
1489 Vc::double_v cmd;
1490
1491 Vc::double_v xm;
1492 Vc::double_v xp;
1493 Vc::double_v ym;
1494 Vc::double_v yp;
1495 Vc::double_v zm;
1496 Vc::double_v zp;
1497
1498 load_crs<x,-1,0>(xm,grid,ids);
1499 load_crs<x,1,0>(xp,grid,ids);
1500 load_crs<y,-1,0>(ym,grid,ids);
1501 load_crs<y,1,0>(yp,grid,ids);
1502 load_crs<z,-1,0>(zm,grid,ids);
1503 load_crs<z,1,0>(zp,grid,ids);
1504 load_crs<x,0,0>(cmd,grid,ids);
1505
1506 Vc::double_v Lap = xm + xp +
1507 ym + yp +
1508 zm + zp - 6.0*cmd;
1509
1510 Vc::Mask<double> surround;
1511
1512 for (int i = 0 ; i < Vc::double_v::Size ; i++)
1513 {surround[i] = (mask_sum[i] == 6);}
1514
1515 Lap = Vc::iif(surround,Lap,Vc::double_v(1.0));
1516
1517 store_crs<1>(grid,Lap,ids);
1518 });
1519 }
1520
1521 int tot_six = 0;
1522 int tot_one = 0;
1523
1524 bool check = true;
1525 auto it2 = grid.getIterator(start,stop);
1526 while (it2.isNext())
1527 {
1528 auto p = it2.get();
1529
1530 check &= (grid.template get<1>(p) == 6 || grid.template get<1>(p) == 1);
1531
1532 // Check the six should be a six
1533 auto xp = p.move(0,1);
1534 auto xm = p.move(0,-1);
1535
1536 auto yp = p.move(1,1);
1537 auto ym = p.move(1,-1);
1538
1539 auto zp = p.move(2,1);
1540 auto zm = p.move(2,-1);
1541
1542 bool is_six;
1543 if (grid.existPoint(xp) && grid.existPoint(xm) &&
1544 grid.existPoint(yp) && grid.existPoint(ym) &&
1545 grid.existPoint(zp) && grid.existPoint(zm))
1546 {
1547 is_six = true;
1548 }
1549 else
1550 {
1551 is_six = false;
1552 }
1553
1554 if (is_six == true && grid.template get<1>(p) != 6.0)
1555 {
1556 check = false;
1557 break;
1558 }
1559
1560 if (grid.template get<1>(p) == 1)
1561 {
1562 tot_one++;
1563 }
1564
1565 if (grid.template get<1>(p) == 6)
1566 {
1567 tot_six++;
1568 }
1569
1570 ++it2;
1571 }
1572
1573 BOOST_REQUIRE_EQUAL(check,true);
1574 BOOST_REQUIRE_EQUAL(tot_six,15857813);
1575 BOOST_REQUIRE_EQUAL(tot_one,2977262);
1576 // Check correct-ness
1577
1578// print_grid("debug_out",grid);
1579}
1580
1581BOOST_AUTO_TEST_CASE( sparse_grid_fast_stencil_vectorized_cross_simplified_ids_vector)
1582{
1583 size_t sz[3] = {501,501,501};
1584 size_t sz_cell[3] = {500,500,500};
1585
1586 sgrid_soa<3,aggregate<double,double,int,double[3],double[3]>,HeapMemory> grid(sz);
1587
1588 grid.getBackgroundValue().template get<0>() = 0.0;
1589
1590 CellDecomposer_sm<3, float, shift<3,float>> cdsm;
1591
1592 Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
1593
1594 cdsm.setDimensions(domain, sz_cell, 0);
1595
1596 fill_sphere_quad_v(grid,cdsm);
1597
1598 //grid.reorder();
1599
1600 grid_key_dx<3> start({1,1,1});
1601 grid_key_dx<3> stop({499,499,499});
1602
1603 for (int i = 0 ; i < 1 ; i++)
1604 {
1605 grid.conv_cross_ids<1,double>(start,stop,[](auto & grid, auto & ids,
1606 unsigned char * mask_sum){
1607 Vc::double_v cmd;
1608
1609 Vc::double_v xm;
1610 Vc::double_v xp;
1611 Vc::double_v ym;
1612 Vc::double_v yp;
1613 Vc::double_v zm;
1614 Vc::double_v zp;
1615
1616 load_crs<x,-1,0>(xm,grid,ids);
1617 load_crs<x,1,0>(xp,grid,ids);
1618 load_crs<y,-1,0>(ym,grid,ids);
1619 load_crs<y,1,0>(yp,grid,ids);
1620 load_crs<z,-1,0>(zm,grid,ids);
1621 load_crs<z,1,0>(zp,grid,ids);
1622 load_crs<x,0,0>(cmd,grid,ids);
1623
1624 Vc::double_v Lap = xm + xp +
1625 ym + yp +
1626 zm + zp - 6.0*cmd;
1627
1628 // Lap for the vector
1629
1630 load_crs_v<x,-1,x,3>(xm,grid,ids);
1631 load_crs_v<x,1,x,3>(xp,grid,ids);
1632 load_crs_v<y,-1,x,3>(ym,grid,ids);
1633 load_crs_v<y,1,x,3>(yp,grid,ids);
1634 load_crs_v<z,-1,x,3>(zm,grid,ids);
1635 load_crs_v<z,1,x,3>(zp,grid,ids);
1636 load_crs_v<x,0,x,3>(cmd,grid,ids);
1637
1638 Vc::double_v Lap_x = xm + xp +
1639 ym + yp +
1640 zm + zp - 6.0*cmd;
1641
1642 load_crs_v<x,-1,y,3>(xm,grid,ids);
1643 load_crs_v<x,1,y,3>(xp,grid,ids);
1644 load_crs_v<y,-1,y,3>(ym,grid,ids);
1645 load_crs_v<y,1,y,3>(yp,grid,ids);
1646 load_crs_v<z,-1,y,3>(zm,grid,ids);
1647 load_crs_v<z,1,y,3>(zp,grid,ids);
1648 load_crs_v<x,0,y,3>(cmd,grid,ids);
1649
1650 Vc::double_v Lap_y = xm + xp +
1651 ym + yp +
1652 zm + zp - 6.0*cmd;
1653
1654 load_crs_v<x,-1,z,3>(xm,grid,ids);
1655 load_crs_v<x,1,z,3>(xp,grid,ids);
1656 load_crs_v<y,-1,z,3>(ym,grid,ids);
1657 load_crs_v<y,1,z,3>(yp,grid,ids);
1658 load_crs_v<z,-1,z,3>(zm,grid,ids);
1659 load_crs_v<z,1,z,3>(zp,grid,ids);
1660 load_crs_v<x,0,z,3>(cmd,grid,ids);
1661
1662 Vc::double_v Lap_z = xm + xp +
1663 ym + yp +
1664 zm + zp - 6.0*cmd;
1665
1666 Vc::Mask<double> surround;
1667
1668 for (int i = 0 ; i < Vc::double_v::Size ; i++)
1669 {surround[i] = (mask_sum[i] == 6);}
1670
1671 Lap = Vc::iif(surround,Lap,Vc::double_v(1.0));
1672 Lap_x = Vc::iif(surround,Lap_x,Vc::double_v(1.0));
1673 Lap_y = Vc::iif(surround,Lap_y,Vc::double_v(1.0));
1674 Lap_z = Vc::iif(surround,Lap_z,Vc::double_v(1.0));
1675
1676 store_crs<1>(grid,Lap,ids);
1677 store_crs_v<4,x>(grid,Lap_x,ids);
1678 store_crs_v<4,y>(grid,Lap_y,ids);
1679 store_crs_v<4,z>(grid,Lap_z,ids);
1680 });
1681 }
1682
1683 int tot_six = 0;
1684 int tot_one = 0;
1685
1686 bool check = true;
1687 auto it2 = grid.getIterator(start,stop);
1688 while (it2.isNext())
1689 {
1690 auto p = it2.get();
1691
1692 check &= (grid.template get<1>(p) == 6 || grid.template get<1>(p) == 1);
1693 check &= (grid.template get<4>(p)[0] == 6 || grid.template get<4>(p)[0] == 1);
1694 check &= (grid.template get<4>(p)[1] == 6 || grid.template get<4>(p)[1] == 1);
1695 check &= (grid.template get<4>(p)[2] == 6 || grid.template get<4>(p)[2] == 1);
1696
1697 // Check the six should be a six
1698 auto xp = p.move(0,1);
1699 auto xm = p.move(0,-1);
1700
1701 auto yp = p.move(1,1);
1702 auto ym = p.move(1,-1);
1703
1704 auto zp = p.move(2,1);
1705 auto zm = p.move(2,-1);
1706
1707 bool is_six;
1708 if (grid.existPoint(xp) && grid.existPoint(xm) &&
1709 grid.existPoint(yp) && grid.existPoint(ym) &&
1710 grid.existPoint(zp) && grid.existPoint(zm))
1711 {
1712 is_six = true;
1713 }
1714 else
1715 {
1716 is_six = false;
1717 }
1718
1719 if (is_six == true && grid.template get<1>(p) != 6.0)
1720 {
1721 check = false;
1722 break;
1723 }
1724
1725 if (is_six == true && grid.template get<4>(p)[0] != 6.0)
1726 {
1727 check = false;
1728 break;
1729 }
1730
1731 if (is_six == true && grid.template get<4>(p)[1] != 6.0)
1732 {
1733 check = false;
1734 break;
1735 }
1736
1737 if (is_six == true && grid.template get<4>(p)[2] != 6.0)
1738 {
1739 check = false;
1740 break;
1741 }
1742
1743 if (grid.template get<1>(p) == 1)
1744 {
1745 tot_one++;
1746 }
1747
1748 if (grid.template get<1>(p) == 6)
1749 {
1750 tot_six++;
1751 }
1752
1753 ++it2;
1754 }
1755
1756 BOOST_REQUIRE_EQUAL(check,true);
1757 BOOST_REQUIRE_EQUAL(tot_six,15857813);
1758 BOOST_REQUIRE_EQUAL(tot_one,2977262);
1759 // Check correct-ness
1760
1761// print_grid("debug_out",grid);
1762}
1763
1764BOOST_AUTO_TEST_CASE( sparse_grid_slow_stencil)
1765{
1766 size_t sz[3] = {501,501,501};
1767 size_t sz_cell[3] = {500,500,500};
1768
1769 sgrid_cpu<3,aggregate<double,double,int>,HeapMemory> grid(sz);
1770
1771 grid.getBackgroundValue().template get<0>() = 0.0;
1772
1773 CellDecomposer_sm<3, float, shift<3,float>> cdsm;
1774
1775 Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
1776
1777 cdsm.setDimensions(domain, sz_cell, 0);
1778
1779 fill_sphere_quad(grid,cdsm);
1780
1781 grid_key_dx<3> start({1,1,1});
1782 grid_key_dx<3> stop({499,499,499});
1783
1784 for (int i = 0 ; i < 1 ; i++)
1785 {
1786
1787 auto it = grid.getIterator();
1788
1789 while (it.isNext())
1790 {
1791 // center point
1792 auto p = it.get();
1793
1794 // plus,minus X,Y,Z
1795 auto mx = p.move(0,-1);
1796 auto px = p.move(0,1);
1797 auto my = p.move(1,-1);
1798 auto py = p.move(1,1);
1799 auto mz = p.move(2,-1);
1800 auto pz = p.move(2,1);
1801
1802 bool surround = grid.existPoint(mx) & grid.existPoint(px) & grid.existPoint(py) & grid.existPoint(my) & grid.existPoint(mz) & grid.existPoint(pz);
1803
1804 double Lap = grid.template get<0>(mz) + grid.template get<0>(pz) +
1805 grid.template get<0>(my) + grid.template get<0>(py) +
1806 grid.template get<0>(mx) + grid.template get<0>(px) - 6.0*grid.template get<0>(p);
1807
1808 grid.template insert<1>(p) = (surround)?Lap:6.0;
1809
1810 ++it;
1811 }
1812
1813 }
1814
1815 bool check = true;
1816 auto it2 = grid.getIterator();
1817 while (it2.isNext())
1818 {
1819 auto p = it2.get();
1820
1821 check &= grid.template get<1>(p) == 6;
1822
1823 ++it2;
1824 }
1825
1826 BOOST_REQUIRE_EQUAL(check,true);
1827 // Check correct-ness
1828 //grid.write("debug_out");
1829}
1830
1831template<typename sgrid> void Test_unpack_and_check_full(sgrid & grid)
1832{
1833 grid_key_dx<3> end;
1834 grid_key_dx<3> zero({0,0,0});
1835 size_t req2 = 0;
1836 size_t req3 = 0;
1837 size_t sz[3];
1838
1839 for (size_t i = 0 ; i < 3 ; i++)
1840 {
1841 end.set_d(i,grid.getGrid().size(i) - 1);
1842 sz[i] = 0;
1843 }
1844
1845 grid.template packRequest<0>(req2);
1846 auto sub_it = grid.getIterator(zero,end);
1847 grid.template packRequest<0>(sub_it,req3);
1848
1849 BOOST_REQUIRE_EQUAL(req2,req3);
1850
1851 Pack_stat sts2;
1852 Pack_stat sts3;
1853
1854 // allocate the memory 2
1855 HeapMemory pmem2;
1856 pmem2.allocate(req2);
1857 ExtPreAlloc<HeapMemory> & mem2 = *(new ExtPreAlloc<HeapMemory>(req2,pmem2));
1858 mem2.incRef();
1859
1860 // allocate the memory 3
1861 HeapMemory pmem3;
1862 pmem3.allocate(req3);
1863 ExtPreAlloc<HeapMemory> & mem3 = *(new ExtPreAlloc<HeapMemory>(req3,pmem3));
1864 mem3.incRef();
1865
1866 grid.template pack<0>(mem2,sts2);
1867 grid.template pack<0>(mem3,sub_it,sts3);
1868
1869 BOOST_REQUIRE_EQUAL(mem2.size(),mem3.size());
1870
1871 bool check = true;
1872
1873 char * p2 = (char *)pmem2.getPointer();
1874 char * p3 = (char *)pmem3.getPointer();
1875
1876 for (size_t i = 0 ; i < mem2.size(); i++)
1877 {
1878 check &= (p2[i] == p3[i]);
1879 }
1880
1881 BOOST_REQUIRE_EQUAL(check,true);
1882
1883 // Unpack on a Sparse grid without information
1884
1885 Unpack_stat ps;
1886 sgrid empty(sz);
1887
1888 empty.template unpack<0>(mem3,ps);
1889
1890 BOOST_REQUIRE_EQUAL(empty.getGrid().size(0),grid.getGrid().size(0));
1891 BOOST_REQUIRE_EQUAL(empty.getGrid().size(1),grid.getGrid().size(1));
1892 BOOST_REQUIRE_EQUAL(empty.getGrid().size(2),grid.getGrid().size(2));
1893
1894 auto it = empty.getIterator();
1895
1896 while (it.isNext())
1897 {
1898 auto p = it.get();
1899
1900 check &= (grid.template get<0>(p) == empty.template get<0>(p));
1901
1902 ++it;
1903 }
1904
1905 BOOST_REQUIRE_EQUAL(check,true);
1906}
1907
1908
1909template<typename sgrid> void Test_unpack_and_check_full_noprp(sgrid & grid)
1910{
1911 grid_key_dx<3> end;
1912 grid_key_dx<3> zero({0,0,0});
1913 size_t req2 = 0;
1914 size_t req3 = 0;
1915 size_t sz[3];
1916
1917 for (size_t i = 0 ; i < 3 ; i++)
1918 {
1919 end.set_d(i,grid.getGrid().size(i) - 1);
1920 sz[i] = 0;
1921 }
1922
1923 grid.template packRequest(req2);
1924 auto sub_it = grid.getIterator(zero,end);
1925 grid.template packRequest<0,1>(sub_it,req3);
1926
1927 BOOST_REQUIRE_EQUAL(req2,req3);
1928
1929 Pack_stat sts2;
1930 Pack_stat sts3;
1931
1932 // allocate the memory 2
1933 HeapMemory pmem2;
1934 pmem2.allocate(req2);
1935 ExtPreAlloc<HeapMemory> & mem2 = *(new ExtPreAlloc<HeapMemory>(req2,pmem2));
1936 mem2.incRef();
1937
1938 // allocate the memory 3
1939 HeapMemory pmem3;
1940 pmem3.allocate(req3);
1941 ExtPreAlloc<HeapMemory> & mem3 = *(new ExtPreAlloc<HeapMemory>(req3,pmem3));
1942 mem3.incRef();
1943
1944 mem2.fill(0);
1945 mem3.fill(0);
1946
1947 grid.template pack(mem2,sts2);
1948 grid.template pack<0,1>(mem3,sub_it,sts3);
1949
1950 BOOST_REQUIRE_EQUAL(mem2.size(),mem3.size());
1951
1952 bool check = true;
1953
1954 char * p2 = (char *)pmem2.getPointer();
1955 char * p3 = (char *)pmem3.getPointer();
1956
1957 for (size_t i = 0 ; i < mem2.size(); i++)
1958 {
1959 check &= (p2[i] == p3[i]);
1960 }
1961
1962 BOOST_REQUIRE_EQUAL(check,true);
1963
1964 // Unpack on a Sparse grid without information
1965
1966 Unpack_stat ps;
1967 sgrid empty(sz);
1968
1969 empty.template unpack(mem2,ps);
1970
1971 BOOST_REQUIRE_EQUAL(empty.getGrid().size(0),grid.getGrid().size(0));
1972 BOOST_REQUIRE_EQUAL(empty.getGrid().size(1),grid.getGrid().size(1));
1973 BOOST_REQUIRE_EQUAL(empty.getGrid().size(2),grid.getGrid().size(2));
1974
1975 auto it = empty.getIterator();
1976
1977 while (it.isNext())
1978 {
1979 auto p = it.get();
1980
1981 check &= (grid.template get<0>(p) == empty.template get<0>(p));
1982
1983 ++it;
1984 }
1985
1986 BOOST_REQUIRE_EQUAL(check,true);
1987}
1988
1989template<typename sgrid> void Test_unpack_and_check(sgrid & grid, sgrid & grid2, size_t (& sz_cell)[3])
1990{
1991 grid.getBackgroundValue().template get<0>() = 0.0;
1992
1993 grid_key_dx<3> start({0,0,0});
1994 grid_key_dx<3> stop({250,250,250});
1995
1996 Box<3,size_t> bx({0,0,0},{250,250,250});
1997
1998 auto sub_it = grid.getIterator(start,stop);
1999 size_t req = 0;
2000 grid.template packRequest<0>(sub_it,req);
2001
2002 // allocate the memory
2003 HeapMemory pmem;
2004 pmem.allocate(req);
2005 ExtPreAlloc<HeapMemory> & mem = *(new ExtPreAlloc<HeapMemory>(req,pmem));
2006 mem.incRef();
2007
2008 Pack_stat sts;
2009
2010 grid.template pack<0>(mem,sub_it,sts);
2011
2012 // now we unpack on another grid
2013
2014 int ctx = 0;
2015 Unpack_stat usts;
2016 grid2.template unpack<0>(mem,sub_it,usts,ctx,rem_copy_opt::NONE_OPT);
2017
2018 bool match = true;
2019 auto it = grid.getIterator();
2020
2021 while (it.isNext())
2022 {
2023 auto key = it.get();
2024
2025 if (bx.isInside(key.toPoint()) == true)
2026 {match &= grid.template get<0>(key) == grid2.template get<0>(key);}
2027
2028 ++it;
2029 }
2030
2031 BOOST_REQUIRE_EQUAL(match,true);
2032}
2033
2034BOOST_AUTO_TEST_CASE( sparse_grid_sub_grid_it_packing)
2035{
2036 size_t sz[3] = {501,501,501};
2037 size_t sz_cell[3] = {500,500,500};
2038
2039 sgrid_cpu<3,aggregate<double,int>,HeapMemory> grid(sz);
2040 sgrid_cpu<3,aggregate<double,int>,HeapMemory> grid2(sz);
2041 sgrid_cpu<3,aggregate<double,int>,HeapMemory> grid3(sz);
2042
2043 CellDecomposer_sm<3, float, shift<3,float>> cdsm;
2044
2045 Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
2046
2047 cdsm.setDimensions(domain, sz_cell, 0);
2048
2049 fill_sphere(grid,cdsm);
2050
2051 Test_unpack_and_check(grid,grid2,sz_cell);
2052
2053 grid_key_dx<3> start({251,251,251});
2054 grid_key_dx<3> stop({490,490,490});
2055
2056
2057 size_t cnt2 = 0;
2058
2059 auto sub_it = grid.getIterator(start,stop);
2060
2061 while (sub_it.isNext())
2062 {
2063 auto p = sub_it.get();
2064
2065 grid3.template insert<0>(p) = grid.template get<0>(p);
2066 grid3.template insert<1>(p) = grid.template get<1>(p);
2067
2068 cnt2++;
2069
2070 ++sub_it;
2071 }
2072
2073 Test_unpack_and_check(grid,grid3,sz_cell);
2074
2075 grid_key_dx<3> start2({251,251,251});
2076 grid_key_dx<3> stop2({490,490,490});
2077
2078 size_t cnt = 0;
2079
2080 auto sub_it2 = grid.getIterator(start2,stop2);
2081
2082 bool match = true;
2083
2084 while (sub_it2.isNext())
2085 {
2086 auto p = sub_it2.get();
2087
2088 match &= grid3.template insert<0>(p) == grid.template get<0>(p);
2089 match &= grid3.template insert<1>(p) == grid.template get<1>(p);
2090
2091 cnt++;
2092
2093 ++sub_it2;
2094 }
2095
2096 BOOST_REQUIRE_EQUAL(match,true);
2097 BOOST_REQUIRE_EQUAL(cnt,cnt2);
2098}
2099
2100
2101BOOST_AUTO_TEST_CASE( sparse_grid_remove_area)
2102{
2103 size_t sz[3] = {501,501,501};
2104
2105 sgrid_cpu<3,aggregate<double,int>,HeapMemory> grid(sz);
2106
2107 Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
2108
2109
2110 Box<3,size_t> bx_create({100,100,100},{400,400,400});
2111 Box<3,long int> bx_delete({150,150,150},{350,350,350});
2112
2113 grid_sm<3,void> gs(sz);
2114 grid_key_dx_iterator_sub<3> sub(gs,bx_create.getKP1(),bx_create.getKP2());
2115
2116 while (sub.isNext())
2117 {
2118 auto p = sub.get();
2119
2120 grid.template insert<0>(p) = 1.0;
2121
2122 ++sub;
2123 }
2124
2125 grid.remove(bx_delete);
2126
2127 bool check = true;
2128
2129 size_t cnt = 0;
2130 auto it2 = grid.getIterator(bx_create.getKP1(),bx_create.getKP2());
2131 while (it2.isNext())
2132 {
2133 auto p = it2.get();
2134
2135 check &= bx_delete.isInside(p.toPoint()) == false;
2136
2137 cnt++;
2138
2139 ++it2;
2140 }
2141
2142 BOOST_REQUIRE_EQUAL(check,true);
2143
2144 BOOST_REQUIRE_EQUAL(cnt,bx_create.getVolumeKey() - bx_delete.getVolumeKey());
2145}
2146
2147BOOST_AUTO_TEST_CASE( sparse_grid_copy_to)
2148{
2149 size_t sz[3] = {501,501,501};
2150
2151 sgrid_cpu<3,aggregate<double,int>,HeapMemory> grid(sz);
2152
2153 grid.getBackgroundValue().template get<0>() = 0.0;
2154
2155 size_t sz_g2[3] = {259,27,27};
2156 sgrid_cpu<3,aggregate<double,int>,HeapMemory> grid2(sz_g2);
2157
2158 grid_key_dx_iterator<3> key_it(grid2.getGrid());
2159 auto gs = grid2.getGrid();
2160
2161 while (key_it.isNext())
2162 {
2163 auto key = key_it.get();
2164
2165 grid2.insert<0>(key) = gs.LinId(key);
2166
2167 ++key_it;
2168 }
2169
2170 Box<3,size_t> bx_src({1,1,1},{255,23,23});
2171 Box<3,size_t> bx_dst({5,5,5},{259,27,27});
2172
2173 grid.copy_to(grid2,bx_src,bx_dst);
2174
2175 BOOST_REQUIRE(grid.size() != 0);
2176 BOOST_REQUIRE_EQUAL(grid.size(),bx_dst.getVolumeKey());
2177
2178 bool match = true;
2179 auto it_check = grid.getIterator(bx_dst.getKP1(),bx_dst.getKP2());
2180
2181 while (it_check.isNext())
2182 {
2183 auto key = it_check.get();
2184
2185 grid_key_dx<3> key2 = key - bx_dst.getKP1();
2186 key2 += bx_src.getKP1();
2187
2188 if (grid.template get<0>(key) != gs.LinId(key2))
2189 {match = false;}
2190
2191 ++it_check;
2192 }
2193
2194 BOOST_REQUIRE_EQUAL(match,true);
2195
2196 bx_dst += Point<3,size_t>({0,40,40});
2197 grid.copy_to(grid2,bx_src,bx_dst);
2198
2199 BOOST_REQUIRE(grid.size() != 0);
2200 BOOST_REQUIRE_EQUAL(grid.size(),2*bx_dst.getVolumeKey());
2201
2202 match = true;
2203 auto it_check2 = grid.getIterator(bx_dst.getKP1(),bx_dst.getKP2());
2204
2205 while (it_check2.isNext())
2206 {
2207 auto key = it_check2.get();
2208
2209 grid_key_dx<3> key2 = key - bx_dst.getKP1();
2210 key2 += bx_src.getKP1();
2211
2212 if (grid.template get<0>(key) != gs.LinId(key2))
2213 {match = false;}
2214
2215 ++it_check2;
2216 }
2217
2218 BOOST_REQUIRE_EQUAL(match,true);
2219}
2220
2221BOOST_AUTO_TEST_CASE( sparse_pack_full )
2222{
2223 size_t sz[3] = {501,501,501};
2224 size_t sz_cell[3] = {500,500,500};
2225
2226 sgrid_cpu<3,aggregate<double,int>,HeapMemory> grid(sz);
2227
2228 CellDecomposer_sm<3, float, shift<3,float>> cdsm;
2229
2230 Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
2231
2232 cdsm.setDimensions(domain, sz_cell, 0);
2233
2234 fill_sphere(grid,cdsm);
2235
2236 Test_unpack_and_check_full(grid);
2237}
2238
2239BOOST_AUTO_TEST_CASE( sparse_pack_full_noprp )
2240{
2241 size_t sz[3] = {501,501,501};
2242 size_t sz_cell[3] = {500,500,500};
2243
2244 sgrid_cpu<3,aggregate<double,int>,HeapMemory> grid(sz);
2245
2246 CellDecomposer_sm<3, float, shift<3,float>> cdsm;
2247
2248 Box<3,float> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
2249
2250 cdsm.setDimensions(domain, sz_cell, 0);
2251
2252 fill_sphere(grid,cdsm);
2253
2254 Test_unpack_and_check_full_noprp(grid);
2255}
2256
2257BOOST_AUTO_TEST_CASE( sparse_operator_equal )
2258{
2259 size_t sz[3] = {270,270,270};
2260
2261 sgrid_cpu<3,aggregate<double,int>,HeapMemory> grid(sz);
2262
2263 grid.getBackgroundValue().template get<0>() = 0.0;
2264
2265 sgrid_cpu<3,aggregate<double,int>,HeapMemory> grid2;
2266
2267 grid_key_dx_iterator<3> key_it(grid.getGrid());
2268 auto gs = grid.getGrid();
2269
2270 while (key_it.isNext())
2271 {
2272 auto key = key_it.get();
2273
2274 grid.insert<0>(key) = gs.LinId(key);
2275
2276 ++key_it;
2277 }
2278
2279 grid2 = grid;
2280
2281 BOOST_REQUIRE(grid.size() == grid2.size());
2282
2283 bool match = true;
2284 auto it_check = grid.getIterator();
2285
2286 while (it_check.isNext())
2287 {
2288 auto key = it_check.get();
2289
2290 if (grid.template get<0>(key) != grid2.template get<0>(key))
2291 {match = false;}
2292
2293 ++it_check;
2294 }
2295
2296 BOOST_REQUIRE_EQUAL(match,true);
2297}
2298
2299BOOST_AUTO_TEST_CASE( sparse_testing_clear )
2300{
2301 size_t sz[3] = {10000,10000,10000};
2302
2303 sgrid_cpu<3,aggregate<float>,HeapMemory> grid(sz);
2304
2305 grid.getBackgroundValue().template get<0>() = 555.0;
2306
2307 grid_key_dx<3> key1({5000,5000,5000});
2308 grid_key_dx<3> key2({5001,5001,5001});
2309 grid_key_dx<3> key3({5002,5003,5003});
2310
2311 grid.template insert<0>(key1) = 1.0;
2312 grid.template insert<0>(key2) = 2.0;
2313 grid.template insert<0>(key3) = 3.0;
2314
2315
2316 grid.clear();
2317
2318 grid_key_dx<3> keyzero({0,0,0});
2319 BOOST_REQUIRE_EQUAL(grid.template get<0>(keyzero),555.0);
2320
2321 grid.template insert<0>(key1) = 1.0;
2322 grid.template insert<0>(key2) = 2.0;
2323 grid.template insert<0>(key3) = 3.0;
2324
2325 BOOST_REQUIRE_EQUAL(grid.template get<0>(key1),1.0);
2326 BOOST_REQUIRE_EQUAL(grid.template get<0>(key2),2.0);
2327 BOOST_REQUIRE_EQUAL(grid.template get<0>(key3),3.0);
2328
2329 BOOST_REQUIRE_EQUAL(grid.template get<0>(keyzero),555.0);
2330}
2331
2332BOOST_AUTO_TEST_SUITE_END()
2333
2334