OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
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 
15 const int x = 0;
16 const int y = 1;
17 const int z = 2;
18 
19 BOOST_AUTO_TEST_SUITE( sgrid_dist_id_test )
20 
21 
22 
23 BOOST_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 
126 BOOST_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 
233 BOOST_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 
340 BOOST_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 
517 BOOST_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 
692 BOOST_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 
868 BOOST_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 
938 BOOST_AUTO_TEST_SUITE_END()
grid_dist_id::write
bool write(std::string output, size_t opt=VTK_WRITER|FORMAT_BINARY)
Write the distributed grid information.
Definition: grid_dist_id.hpp:2886
grid_key_dx
grid_key_dx is the key to access any element in the grid
Definition: grid_key.hpp:18
grid_dist_id::get
auto get(const grid_dist_key_dx< dim, bg_key > &v1) const -> typename std::add_lvalue_reference< decltype(loc_grid.get(v1.getSub()).template get< p >(v1.getKey()))>::type
Get the reference of the selected element.
Definition: grid_dist_id.hpp:2221
grid_dist_key_dx
Grid key for a distributed grid.
Definition: grid_dist_key.hpp:6
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
Ghost
Definition: Ghost.hpp:39
grid_dist_id
This is a distributed grid.
Definition: grid_dist_id.hpp:96
cross_stencil_v
Definition: SparseGrid_conv_opt.hpp:204
grid_dist_id::getGridIterator
grid_dist_id_iterator_dec< Decomposition > getGridIterator(const grid_key_dx< dim > &start, const grid_key_dx< dim > &stop)
Definition: grid_dist_id.hpp:1787
cub::int
KeyT const ValueT ValueT OffsetIteratorT OffsetIteratorT int
[in] The number of segments that comprise the sorting data
Definition: dispatch_radix_sort.cuh:331
Box
This class represent an N-dimensional box.
Definition: Box.hpp:60
periodicity
Boundary conditions.
Definition: common.hpp:21
F
[v_transform metafunction]
Definition: variadic_to_vmpl_unit_test.hpp:17