OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
sgrid_dist_id_unit_tests.cpp
1/*
2 * sgrid_dist_id_unit_tests.cpp
3 *
4 * Created on: Nov 18, 2017
5 * Author: i-bird
6 */
7
8
9#define BOOST_TEST_DYN_LINK
10#include <boost/test/unit_test.hpp>
11#include "Grid/grid_dist_id.hpp"
12#include "Point_test.hpp"
13
14
15const int x = 0;
16const int y = 1;
17const int z = 2;
18
19BOOST_AUTO_TEST_SUITE( sgrid_dist_id_test )
20
21
22
23BOOST_AUTO_TEST_CASE (sgrid_dist_id_soa )
24{
25 periodicity<3> bc = {PERIODIC, PERIODIC, PERIODIC};
26
27 // Domain
28 Box<3,double> domain({-0.3,-0.3,-0.3},{1.0,1.0,1.0});
29
30 // grid size
31 size_t sz[3];
32 sz[0] = 1024;
33 sz[1] = 1024;
34 sz[2] = 1024;
35
36 // Ghost
37 Ghost<3,double> g(0.01);
38
39 sgrid_dist_soa<3,double,Point_test<float>> sg(sz,domain,g,bc);
40
41 // create a grid iterator over a bilion point
42
43 auto it = sg.getGridIterator();
44
45 while(it.isNext())
46 {
47 auto gkey = it.get();
48 auto key = it.get_dist();
49
50 size_t sx = gkey.get(0) - 512;
51 size_t sy = gkey.get(1) - 512;
52 size_t sz = gkey.get(2) - 512;
53
54 if (sx*sx + sy*sy + sz*sz < 128*128)
55 {
56 sg.template insert<0>(key) = 1.0;
57 }
58
59 ++it;
60 }
61
62 bool match = true;
63 auto it2 = sg.getGridIterator();
64
65 while(it2.isNext())
66 {
67 auto gkey = it2.get();
68 auto key = it2.get_dist();
69
70 size_t sx = gkey.get(0) - 512;
71 size_t sy = gkey.get(1) - 512;
72 size_t sz = gkey.get(2) - 512;
73
74 if (sx*sx + sy*sy + sz*sz < 128*128)
75 {
76 match &= (sg.template get<0>(key) == 1.0);
77 }
78
79 ++it2;
80 }
81
82 auto & gr = sg.getGridInfo();
83
84 auto it3 = sg.getDomainIterator();
85
86 while (it3.isNext())
87 {
88 auto key = it3.get();
89 auto gkey = it3.getGKey(key);
90
91 sg.template insert<0>(key) = gkey.get(0)*gkey.get(0) + gkey.get(1)*gkey.get(1) + gkey.get(2)*gkey.get(2);
92
93 ++it3;
94 }
95
96 sg.ghost_get<0>();
97 // now we check the stencil
98
99 bool good = true;
100 auto it4 = sg.getDomainIterator();
101
102 while (it4.isNext())
103 {
104 auto key = it4.get();
105 auto gkey = it4.getGKey(key);
106
107 size_t sx = gkey.get(0) - 512;
108 size_t sy = gkey.get(1) - 512;
109 size_t sz = gkey.get(2) - 512;
110
111 double lap;
112
113 lap = sg.template get<0>(key.move(x,1)) + sg.template get<0>(key.move(x,-1)) +
114 sg.template get<0>(key.move(y,1)) + sg.template get<0>(key.move(y,-1)) +
115 sg.template get<0>(key.move(z,1)) + sg.template get<0>(key.move(z,-1)) -
116 6.0*sg.template get<0>(key);
117
118 good &= (lap == 6.0);
119
120 ++it4;
121 }
122
123 BOOST_REQUIRE_EQUAL(match,true);
124}
125
126BOOST_AUTO_TEST_CASE( sgrid_dist_id_basic_test_2D)
127{
128 periodicity<2> bc = {NON_PERIODIC, NON_PERIODIC};
129
130 // Domain
131 Box<2,double> domain({-0.3,-0.3},{1.0,1.0});
132
133 // grid size
134 size_t sz[2];
135 sz[0] = 1024;
136 sz[1] = 1024;
137
138 // Ghost
139 Ghost<2,double> g(0.01);
140
141 sgrid_dist_id<2,double,Point_test<double>> sg(sz,domain,g,bc);
142
143 // create a grid iterator
144
145 auto it = sg.getGridIterator();
146
147 while(it.isNext())
148 {
149 auto gkey = it.get();
150 auto key = it.get_dist();
151
152
153 long int sx = gkey.get(0) - 512;
154 long int sy = gkey.get(1) - 512;
155
156 if (sx*sx + sy*sy < 128*128)
157 {
158 sg.template insert<0>(key) = 1.0;
159 }
160
161 ++it;
162 }
163
164 bool match = true;
165 auto it2 = sg.getGridIterator();
166
167 while(it2.isNext())
168 {
169 auto gkey = it2.get();
170 auto key = it2.get_dist();
171
172 long int sx = gkey.get(0) - 512;
173 long int sy = gkey.get(1) - 512;
174
175 if (sx*sx + sy*sy < 128*128)
176 {
177 match &= (sg.template get<0>(key) == 1.0);
178 }
179
180 ++it2;
181 }
182
183 BOOST_REQUIRE_EQUAL(match,true);
184
185 auto & gr = sg.getGridInfo();
186
187 auto it3 = sg.getDomainIterator();
188
189 while (it3.isNext())
190 {
191 auto key = it3.get();
192 auto gkey = it3.getGKey(key);
193
194 sg.template insert<0>(key) = gkey.get(0)*gkey.get(0) + gkey.get(1)*gkey.get(1);
195
196 ++it3;
197 }
198
199 sg.ghost_get<0>();
200
201 // now we check the stencil
202
203 bool good = true;
204 auto it4 = sg.getDomainIterator();
205
206 while (it4.isNext())
207 {
208 auto key = it4.get();
209 auto gkey = it4.getGKey(key);
210
211 double lap;
212
213 // Here we check that all point of the stencil are inside*/
214
215 long int sx = gkey.get(0) - 512;
216 long int sy = gkey.get(1) - 512;
217
218 if (sx*sx + sy*sy < 126*126)
219 {
220 lap = sg.template get<0>(key.move(x,1)) + sg.template get<0>(key.move(x,-1)) +
221 sg.template get<0>(key.move(y,1)) + sg.template get<0>(key.move(y,-1)) -
222 4.0*sg.template get<0>(key);
223
224 good &= (lap == 4.0);
225 }
226
227 ++it4;
228 }
229
230 BOOST_REQUIRE_EQUAL(good,true);
231}
232
233BOOST_AUTO_TEST_CASE( sgrid_dist_id_basic_test)
234{
235 periodicity<3> bc = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
236
237 // Domain
238 Box<3,double> domain({-0.3,-0.3,-0.3},{1.0,1.0,1.0});
239
240 // grid size
241 size_t sz[3];
242 sz[0] = 1024;
243 sz[1] = 1024;
244 sz[2] = 1024;
245
246 // Ghost
247 Ghost<3,double> g(0.01);
248
249 sgrid_dist_id<3,double,Point_test<float>> sg(sz,domain,g,bc);
250
251 // create a grid iterator over a bilion point
252
253 auto it = sg.getGridIterator();
254
255 while(it.isNext())
256 {
257 auto gkey = it.get();
258 auto key = it.get_dist();
259
260 size_t sx = gkey.get(0) - 512;
261 size_t sy = gkey.get(1) - 512;
262 size_t sz = gkey.get(2) - 512;
263
264 if (sx*sx + sy*sy + sz*sz < 128*128)
265 {
266 sg.template insert<0>(key) = 1.0;
267 }
268
269 ++it;
270 }
271
272 bool match = true;
273 auto it2 = sg.getGridIterator();
274
275 while(it2.isNext())
276 {
277 auto gkey = it2.get();
278 auto key = it2.get_dist();
279
280 size_t sx = gkey.get(0) - 512;
281 size_t sy = gkey.get(1) - 512;
282 size_t sz = gkey.get(2) - 512;
283
284 if (sx*sx + sy*sy + sz*sz < 128*128)
285 {
286 match &= (sg.template get<0>(key) == 1.0);
287 }
288
289 ++it2;
290 }
291
292 auto & gr = sg.getGridInfo();
293
294 auto it3 = sg.getDomainIterator();
295
296 while (it3.isNext())
297 {
298 auto key = it3.get();
299 auto gkey = it3.getGKey(key);
300
301 sg.template insert<0>(key) = gkey.get(0)*gkey.get(0) + gkey.get(1)*gkey.get(1) + gkey.get(2)*gkey.get(2);
302
303 ++it3;
304 }
305
306 sg.ghost_get<0>();
307 // now we check the stencil
308
309 bool good = true;
310 auto it4 = sg.getDomainIterator();
311
312 while (it4.isNext())
313 {
314 auto key = it4.get();
315 auto gkey = it4.getGKey(key);
316
317 size_t sx = gkey.get(0) - 512;
318 size_t sy = gkey.get(1) - 512;
319 size_t sz = gkey.get(2) - 512;
320
321 if (sx*sx + sy*sy + sz*sz < 125*125)
322 {
323 double lap;
324
325 lap = sg.template get<0>(key.move(x,1)) + sg.template get<0>(key.move(x,-1)) +
326 sg.template get<0>(key.move(y,1)) + sg.template get<0>(key.move(y,-1)) +
327 sg.template get<0>(key.move(z,1)) + sg.template get<0>(key.move(z,-1)) -
328 6.0*sg.template get<0>(key);
329
330 good &= (lap == 6.0);
331 }
332
333 ++it4;
334 }
335
336 BOOST_REQUIRE_EQUAL(match,true);
337}
338
339
340BOOST_AUTO_TEST_CASE( sparse_grid_fast_stencil_vectorized_simplified_conv2)
341{
342 constexpr int U = 0;
343 constexpr int V = 1;
344
345 constexpr int U_next = 2;
346 constexpr int V_next = 3;
347
348 constexpr int x = 0;
349 constexpr int y = 1;
350 constexpr int z = 2;
351
352 Box<3,double> domain({0.0,0.0,0.0},{2.5,2.5,2.5});
353
354 // grid size
355 size_t sz[3] = {32,32,32};
356
357 // Define periodicity of the grid
358 periodicity<3> bc = {PERIODIC,PERIODIC,PERIODIC};
359
360 // Ghost in grid unit
361 Ghost<3,long int> g(1);
362
363 // deltaT
364 double deltaT = 1;
365
366 // Diffusion constant for specie U
367 double du = 2*1e-5;
368
369 // Diffusion constant for specie V
370 double dv = 1*1e-5;
371
372 // Number of timesteps
373 size_t timeSteps = 5000;
374
375 // K and F (Physical constant in the equation)
376 double K = 0.053;
377 double F = 0.014;
378
379 sgrid_dist_soa<3, double, aggregate<double,double,double,double>> grid(sz,domain,g,bc);
380
381 auto it = grid.getGridIterator();
382
383 while (it.isNext())
384 {
385 // Get the local grid key
386 auto key = it.get_dist();
387
388 // Old values U and V
389 grid.template insert<U>(key) = 1.0;
390 grid.template insert<V>(key) = 0.0;
391
392 // Old values U and V
393 grid.template insert<U_next>(key) = 0.0;
394 grid.template insert<V_next>(key) = 0.0;
395
396 ++it;
397 }
398
399 long int x_start = grid.size(0)*1.55f/domain.getHigh(0);
400 long int y_start = grid.size(1)*1.55f/domain.getHigh(1);
401 long int z_start = grid.size(1)*1.55f/domain.getHigh(2);
402
403 long int x_stop = grid.size(0)*1.85f/domain.getHigh(0);
404 long int y_stop = grid.size(1)*1.85f/domain.getHigh(1);
405 long int z_stop = grid.size(1)*1.85f/domain.getHigh(2);
406
407 grid_key_dx<3> start({x_start,y_start,z_start});
408 grid_key_dx<3> stop ({x_stop,y_stop,z_stop});
409 auto it_init = grid.getGridIterator(start,stop);
410
411 while (it_init.isNext())
412 {
413 auto key = it_init.get_dist();
414
415 grid.template insert<U>(key) = 0.5 + (((double)std::rand())/RAND_MAX -0.5)/10.0;
416 grid.template insert<V>(key) = 0.25 + (((double)std::rand())/RAND_MAX -0.5)/20.0;
417
418 ++it_init;
419 }
420
421 // spacing of the grid on x and y
422 double spacing[3] = {grid.spacing(0),grid.spacing(1),grid.spacing(2)};
423 // sync the ghost
424 size_t count = 0;
425 grid.template ghost_get<U,V>();
426
427 // because we assume that spacing[x] == spacing[y] we use formula 2
428 // and we calculate the prefactor of Eq 2
429 double uFactor = deltaT * du/(spacing[x]*spacing[x]);
430 double vFactor = deltaT * dv/(spacing[x]*spacing[x]);
431
432 int stencil[6][3] = {{1,0,0},{-1,0,0},{0,-1,0},{0,1,0},{0,0,-1},{0,0,1}};
433
434
436
437
438 auto func = [uFactor,vFactor,deltaT,F,K](Vc::double_v & u_out,Vc::double_v & v_out,
439 Vc::double_v (& u)[7],Vc::double_v (& v)[7],
440 unsigned char * mask){
441
442 u_out = u[0] + uFactor *(u[1] + u[2] +
443 u[3] + u[4] +
444 u[5] + u[6] - 6.0*u[0]) - deltaT * u[0]*v[0]*v[0]
445 - deltaT * F * (u[0] - 1.0);
446
447 v_out = v[0] + vFactor *(v[1] + v[2] +
448 v[3] + v[4] +
449 v[5] + v[6] - 6.0*v[0]) + deltaT * u[0]*v[0]*v[0]
450 - deltaT * (F+K) * v[0];
451 };
452
453 grid.conv2<U,V,U_next,V_next,1>(stencil,{0,0,0},{(long int)sz[0]-1,(long int)sz[1]-1,(long int)sz[2]-1},func);
454 grid.conv2<U,V,U_next,V_next,1>(stencil,{0,0,0},{(long int)sz[0]-1,(long int)sz[1]-1,(long int)sz[2]-1},func);
455
456 bool match = true;
457
458 {
459 auto it = grid.getDomainIterator();
460
461 double max_U = 0.0;
462 double max_V = 0.0;
463 grid_dist_key_dx<3> k_max;
464 while (it.isNext())
465 {
466 // center point
467 auto Cp = it.get();
468
469 // plus,minus X,Y,Z
470 auto mx = Cp.move(0,-1);
471 auto px = Cp.move(0,+1);
472 auto my = Cp.move(1,-1);
473 auto py = Cp.move(1,1);
474 auto mz = Cp.move(2,-1);
475 auto pz = Cp.move(2,1);
476
477 // update based on Eq 2
478 if ( fabs(grid.get<U>(Cp) + uFactor * (
479 grid.get<U>(mz) +
480 grid.get<U>(pz) +
481 grid.get<U>(my) +
482 grid.get<U>(py) +
483 grid.get<U>(mx) +
484 grid.get<U>(px) -
485 6.0*grid.get<U>(Cp)) +
486 - deltaT * grid.get<U>(Cp) * grid.get<V>(Cp) * grid.get<V>(Cp) +
487 - deltaT * F * (grid.get<U>(Cp) - 1.0) - grid.get<U_next>(Cp)) > 0.000000001 )
488 {
489 match = false;
490 break;
491 }
492
493 // update based on Eq 2
494 if ( fabs(grid.get<V>(Cp) + vFactor * (
495 grid.get<V>(mz) +
496 grid.get<V>(pz) +
497 grid.get<V>(my) +
498 grid.get<V>(py) +
499 grid.get<V>(mx) +
500 grid.get<V>(px) -
501 6*grid.get<V>(Cp)) +
502 deltaT * grid.get<U>(Cp) * grid.get<V>(Cp) * grid.get<V>(Cp) +
503 - deltaT * (F+K) * grid.get<V>(Cp) - grid.get<V_next>(Cp)) > 0.000000001 )
504 {
505 match = false;
506 break;
507 }
508
509 ++it;
510 }
511 }
512
513 BOOST_REQUIRE_EQUAL(match,true);
514}
515
516
517BOOST_AUTO_TEST_CASE( sparse_grid_fast_stencil_vectorized_simplified_conv2_crossing)
518{
519 constexpr int U = 0;
520 constexpr int V = 1;
521
522 constexpr int U_next = 2;
523 constexpr int V_next = 3;
524
525 constexpr int x = 0;
526 constexpr int y = 1;
527 constexpr int z = 2;
528
529 Box<3,double> domain({0.0,0.0,0.0},{2.5,2.5,2.5});
530
531 // grid size
532 size_t sz[3] = {32,32,32};
533
534 // Define periodicity of the grid
535 periodicity<3> bc = {PERIODIC,PERIODIC,PERIODIC};
536
537 // Ghost in grid unit
538 Ghost<3,long int> g(1);
539
540 // deltaT
541 double deltaT = 1;
542
543 // Diffusion constant for specie U
544 double du = 2*1e-5;
545
546 // Diffusion constant for specie V
547 double dv = 1*1e-5;
548
549 // Number of timesteps
550 size_t timeSteps = 5000;
551
552 // K and F (Physical constant in the equation)
553 double K = 0.053;
554 double F = 0.014;
555
556 sgrid_dist_soa<3, double, aggregate<double,double,double,double>> grid(sz,domain,g,bc);
557
558 auto it = grid.getGridIterator();
559
560 while (it.isNext())
561 {
562 // Get the local grid key
563 auto key = it.get_dist();
564
565 // Old values U and V
566 grid.template insert<U>(key) = 1.0;
567 grid.template insert<V>(key) = 0.0;
568
569 // Old values U and V
570 grid.template insert<U_next>(key) = 0.0;
571 grid.template insert<V_next>(key) = 0.0;
572
573 ++it;
574 }
575
576 long int x_start = grid.size(0)*1.55f/domain.getHigh(0);
577 long int y_start = grid.size(1)*1.55f/domain.getHigh(1);
578 long int z_start = grid.size(1)*1.55f/domain.getHigh(2);
579
580 long int x_stop = grid.size(0)*1.85f/domain.getHigh(0);
581 long int y_stop = grid.size(1)*1.85f/domain.getHigh(1);
582 long int z_stop = grid.size(1)*1.85f/domain.getHigh(2);
583
584 grid_key_dx<3> start({x_start,y_start,z_start});
585 grid_key_dx<3> stop ({x_stop,y_stop,z_stop});
586 auto it_init = grid.getGridIterator(start,stop);
587
588 while (it_init.isNext())
589 {
590 auto key = it_init.get_dist();
591
592 grid.template insert<U>(key) = 0.5 + (((double)std::rand())/RAND_MAX -0.5)/10.0;
593 grid.template insert<V>(key) = 0.25 + (((double)std::rand())/RAND_MAX -0.5)/20.0;
594
595 ++it_init;
596 }
597
598 // spacing of the grid on x and y
599 double spacing[3] = {grid.spacing(0),grid.spacing(1),grid.spacing(2)};
600 // sync the ghost
601 size_t count = 0;
602 grid.template ghost_get<U,V>();
603
604 // because we assume that spacing[x] == spacing[y] we use formula 2
605 // and we calculate the prefactor of Eq 2
606 double uFactor = deltaT * du/(spacing[x]*spacing[x]);
607 double vFactor = deltaT * dv/(spacing[x]*spacing[x]);
608
609
611
612
613 auto func = [uFactor,vFactor,deltaT,F,K](Vc::double_v & u_out,Vc::double_v & v_out,
614 Vc::double_v & u,Vc::double_v & v,
615 cross_stencil_v<double> & us,cross_stencil_v<double> & vs,
616 unsigned char * mask){
617
618 u_out = u + uFactor *(us.xm + us.xp +
619 us.ym + us.yp +
620 us.zm + us.zp - 6.0*u) - deltaT * u*v*v
621 - deltaT * F * (u - 1.0);
622
623 v_out = v + vFactor *(vs.xm + vs.xp +
624 vs.ym + vs.yp +
625 vs.zm + vs.zp - 6.0*v) + deltaT * u*v*v
626 - deltaT * (F+K) * v;
627 };
628
629 grid.conv_cross2<U,V,U_next,V_next,1>({0,0,0},{(long int)sz[0]-1,(long int)sz[1]-1,(long int)sz[2]-1},func);
630 grid.conv_cross2<U,V,U_next,V_next,1>({0,0,0},{(long int)sz[0]-1,(long int)sz[1]-1,(long int)sz[2]-1},func);
631
632 bool match = true;
633
634 {
635 auto it = grid.getDomainIterator();
636
637 double max_U = 0.0;
638 double max_V = 0.0;
639 grid_dist_key_dx<3> k_max;
640 while (it.isNext())
641 {
642 // center point
643 auto Cp = it.get();
644
645 // plus,minus X,Y,Z
646 auto mx = Cp.move(0,-1);
647 auto px = Cp.move(0,+1);
648 auto my = Cp.move(1,-1);
649 auto py = Cp.move(1,1);
650 auto mz = Cp.move(2,-1);
651 auto pz = Cp.move(2,1);
652
653 // update based on Eq 2
654 if ( fabs(grid.get<U>(Cp) + uFactor * (
655 grid.get<U>(mz) +
656 grid.get<U>(pz) +
657 grid.get<U>(my) +
658 grid.get<U>(py) +
659 grid.get<U>(mx) +
660 grid.get<U>(px) -
661 6.0*grid.get<U>(Cp)) +
662 - deltaT * grid.get<U>(Cp) * grid.get<V>(Cp) * grid.get<V>(Cp) +
663 - deltaT * F * (grid.get<U>(Cp) - 1.0) - grid.get<U_next>(Cp)) > 0.000000001 )
664 {
665 match = false;
666 break;
667 }
668
669 // update based on Eq 2
670 if ( fabs(grid.get<V>(Cp) + vFactor * (
671 grid.get<V>(mz) +
672 grid.get<V>(pz) +
673 grid.get<V>(my) +
674 grid.get<V>(py) +
675 grid.get<V>(mx) +
676 grid.get<V>(px) -
677 6*grid.get<V>(Cp)) +
678 deltaT * grid.get<U>(Cp) * grid.get<V>(Cp) * grid.get<V>(Cp) +
679 - deltaT * (F+K) * grid.get<V>(Cp) - grid.get<V_next>(Cp)) > 0.000000001 )
680 {
681 match = false;
682 break;
683 }
684
685 ++it;
686 }
687 }
688
689 BOOST_REQUIRE_EQUAL(match,true);
690}
691
692BOOST_AUTO_TEST_CASE( sparse_grid_fast_stencil_vectorized_simplified_conv2_crossing_float)
693{
694 constexpr int U = 0;
695 constexpr int V = 1;
696
697 constexpr int U_next = 2;
698 constexpr int V_next = 3;
699
700 constexpr int x = 0;
701 constexpr int y = 1;
702 constexpr int z = 2;
703
704 Box<3,float> domain({0.0,0.0,0.0},{2.5,2.5,2.5});
705
706 // grid size
707 size_t sz[3] = {32,32,32};
708
709 // Define periodicity of the grid
710 periodicity<3> bc = {PERIODIC,PERIODIC,PERIODIC};
711
712 // Ghost in grid unit
713 Ghost<3,long int> g(1);
714
715 // deltaT
716 float deltaT = 1;
717
718 // Diffusion constant for specie U
719 float du = 2*1e-5;
720
721 // Diffusion constant for specie V
722 float dv = 1*1e-5;
723
724 // Number of timesteps
725 size_t timeSteps = 5000;
726
727 // K and F (Physical constant in the equation)
728 float K = 0.053;
729 float F = 0.014;
730
731 sgrid_dist_soa<3, float, aggregate<float,float,float,float>> grid(sz,domain,g,bc);
732
733 auto it = grid.getGridIterator();
734
735 while (it.isNext())
736 {
737 // Get the local grid key
738 auto key = it.get_dist();
739
740 // Old values U and V
741 grid.template insert<U>(key) = 1.0;
742 grid.template insert<V>(key) = 0.0;
743
744 // Old values U and V
745 grid.template insert<U_next>(key) = 0.0;
746 grid.template insert<V_next>(key) = 0.0;
747
748 ++it;
749 }
750
751 long int x_start = grid.size(0)*1.55f/domain.getHigh(0);
752 long int y_start = grid.size(1)*1.55f/domain.getHigh(1);
753 long int z_start = grid.size(1)*1.55f/domain.getHigh(2);
754
755 long int x_stop = grid.size(0)*1.85f/domain.getHigh(0);
756 long int y_stop = grid.size(1)*1.85f/domain.getHigh(1);
757 long int z_stop = grid.size(1)*1.85f/domain.getHigh(2);
758
759 grid_key_dx<3> start({x_start,y_start,z_start});
760 grid_key_dx<3> stop ({x_stop,y_stop,z_stop});
761 auto it_init = grid.getGridIterator(start,stop);
762
763 while (it_init.isNext())
764 {
765 auto key = it_init.get_dist();
766
767 grid.template insert<U>(key) = 0.5 + (((float)std::rand())/RAND_MAX -0.5)/10.0;
768 grid.template insert<V>(key) = 0.25 + (((float)std::rand())/RAND_MAX -0.5)/20.0;
769
770 ++it_init;
771 }
772
773 // spacing of the grid on x and y
774 float spacing[3] = {grid.spacing(0),grid.spacing(1),grid.spacing(2)};
775 // sync the ghost
776 size_t count = 0;
777 grid.template ghost_get<U,V>();
778
779 // because we assume that spacing[x] == spacing[y] we use formula 2
780 // and we calculate the prefactor of Eq 2
781 float uFactor = deltaT * du/(spacing[x]*spacing[x]);
782 float vFactor = deltaT * dv/(spacing[x]*spacing[x]);
783
784
786
787
788 auto func = [uFactor,vFactor,deltaT,F,K](Vc::float_v & u_out,Vc::float_v & v_out,
789 Vc::float_v & u,Vc::float_v & v,
790 cross_stencil_v<float> & us,cross_stencil_v<float> & vs,
791 unsigned char * mask){
792
793 u_out = u + uFactor *(us.xm + us.xp +
794 us.ym + us.yp +
795 us.zm + us.zp - 6.0f*u) - deltaT * u*v*v
796 - deltaT * F * (u - 1.0f);
797
798 v_out = v + vFactor *(vs.xm + vs.xp +
799 vs.ym + vs.yp +
800 vs.zm + vs.zp - 6.0f*v) + deltaT * u*v*v
801 - deltaT * (F+K) * v;
802 };
803
804 grid.conv_cross2<U,V,U_next,V_next,1>({0,0,0},{(long int)sz[0]-1,(long int)sz[1]-1,(long int)sz[2]-1},func);
805 grid.conv_cross2<U,V,U_next,V_next,1>({0,0,0},{(long int)sz[0]-1,(long int)sz[1]-1,(long int)sz[2]-1},func);
806
807 bool match = true;
808
809 {
810 auto it = grid.getDomainIterator();
811
812 float max_U = 0.0;
813 float max_V = 0.0;
814 grid_dist_key_dx<3> k_max;
815 while (it.isNext())
816 {
817 // center point
818 auto Cp = it.get();
819
820 // plus,minus X,Y,Z
821 auto mx = Cp.move(0,-1);
822 auto px = Cp.move(0,+1);
823 auto my = Cp.move(1,-1);
824 auto py = Cp.move(1,1);
825 auto mz = Cp.move(2,-1);
826 auto pz = Cp.move(2,1);
827
828 // update based on Eq 2
829 if ( fabs(grid.get<U>(Cp) + uFactor * (
830 grid.get<U>(mz) +
831 grid.get<U>(pz) +
832 grid.get<U>(my) +
833 grid.get<U>(py) +
834 grid.get<U>(mx) +
835 grid.get<U>(px) -
836 6.0*grid.get<U>(Cp)) +
837 - deltaT * grid.get<U>(Cp) * grid.get<V>(Cp) * grid.get<V>(Cp) +
838 - deltaT * F * (grid.get<U>(Cp) - 1.0) - grid.get<U_next>(Cp)) > 0.00001 )
839 {
840 match = false;
841 break;
842 }
843
844 // update based on Eq 2
845 if ( fabs(grid.get<V>(Cp) + vFactor * (
846 grid.get<V>(mz) +
847 grid.get<V>(pz) +
848 grid.get<V>(my) +
849 grid.get<V>(py) +
850 grid.get<V>(mx) +
851 grid.get<V>(px) -
852 6*grid.get<V>(Cp)) +
853 deltaT * grid.get<U>(Cp) * grid.get<V>(Cp) * grid.get<V>(Cp) +
854 - deltaT * (F+K) * grid.get<V>(Cp) - grid.get<V_next>(Cp)) > 0.00001 )
855 {
856 match = false;
857 break;
858 }
859
860 ++it;
861 }
862 }
863
864 BOOST_REQUIRE_EQUAL(match,true);
865}
866
867
868BOOST_AUTO_TEST_CASE (sgrid_dist_id_soa_write )
869{
870 periodicity<3> bc = {PERIODIC, PERIODIC, PERIODIC};
871
872 auto & v_cl = create_vcluster<>();
873
874 if (v_cl.size() > 16)
875 {return;}
876
877 // Domain
878 Box<3,double> domain({-0.3,-0.3,-0.3},{1.0,1.0,1.0});
879
880 // grid size
881 size_t sz[3];
882 sz[0] = 256;
883 sz[1] = 256;
884 sz[2] = 256;
885
886 // Ghost
887 Ghost<3,long int> g(1);
888
889 sgrid_dist_soa<3,double,aggregate<double,double[3]>> sg1(sz,domain,g,bc);
890 sgrid_dist_id<3,double,aggregate<double,double[3]>> sg2(sg1.getDecomposition(),sz,g);
891
892 // create a grid iterator over a bilion point
893
894 auto it = sg1.getGridIterator();
895
896 while(it.isNext())
897 {
898 auto gkey = it.get();
899 auto key = it.get_dist();
900
901 size_t sx = gkey.get(0) - 128;
902 size_t sy = gkey.get(1) - 128;
903 size_t sz = gkey.get(2) - 128;
904
905 if (sx*sx + sy*sy + sz*sz < 32*32)
906 {
907 sg1.template insert<0>(key) = 1.0;
908 sg1.template insert<1>(key)[0] = gkey.get(0);
909 sg1.template insert<1>(key)[1] = gkey.get(1);
910 sg1.template insert<1>(key)[2] = gkey.get(2);
911
912 sg2.template insert<0>(key) = 1.0;
913 sg2.template insert<1>(key)[0] = gkey.get(0);
914 sg2.template insert<1>(key)[1] = gkey.get(1);
915 sg2.template insert<1>(key)[2] = gkey.get(2);
916 }
917
918 ++it;
919 }
920
921 sg1.write("sg1_test");
922 sg2.write("sg2_test");
923
924 bool test = compare("sg1_test_" + std::to_string(v_cl.rank()) + ".vtk","sg2_test_" + std::to_string(v_cl.rank()) + ".vtk");
925 BOOST_REQUIRE_EQUAL(true,test);
926
927 sg1.save("hdf5_w1_test");
928 sg2.save("hdf5_w2_test");
929
930 // To uncomment and check
931// sgrid_dist_soa<3,double,aggregate<double,double[3]>> sg1_(sz,domain,g,bc);
932// sgrid_dist_id<3,double,aggregate<double,double[3]>> sg2_(sg1.getDecomposition(),sz,g);
933
934// sg1.load("hdf5_w1_test");
935// sg2.load("hdf5_w2_test");
936}
937
938BOOST_AUTO_TEST_SUITE_END()
Box
This class represent an N-dimensional box.
Definition Box.hpp:61
grid_dist_id
This is a distributed grid.
Definition grid_dist_id.hpp:278
grid_dist_key_dx
Grid key for a distributed grid.
Definition grid_dist_key.hpp:44
grid_dist_key_dx::move
grid_dist_key_dx< dim, base_key > move(size_t i, int s) const
Create a new key moving the old one.
Definition grid_dist_key.hpp:132
grid_key_dx
grid_key_dx is the key to access any element in the grid
Definition grid_key.hpp:19
cub::int
KeyT const ValueT ValueT OffsetIteratorT OffsetIteratorT int
[in] The number of segments that comprise the sorting data
Definition dispatch_radix_sort.cuh:336
F
[v_transform metafunction]
Definition variadic_to_vmpl_unit_test.hpp:18
cross_stencil_v
Definition SparseGrid_conv_opt.hpp:205
periodicity
Boundary conditions.
Definition common.hpp:22