1/*
2 * interpolation_unit_tests.hpp
3 *
4 * Created on: May 5, 2017
5 * Author: i-bird
6 */
7
8#ifndef OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_UNIT_TESTS_HPP_
9#define OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_UNIT_TESTS_HPP_
10
11#define BOOST_TEST_DYN_LINK
12#include <boost/test/unit_test.hpp>
13
14#include "interpolation/mp4_kernel.hpp"
15#include "interpolation/z_spline.hpp"
16#include "interpolation.hpp"
17#include <boost/math/special_functions/pow.hpp>
18#include <Vector/vector_dist.hpp>
19#include <Grid/grid_dist_id.hpp>
20
21BOOST_AUTO_TEST_SUITE( interpolation_test )
22
23template<typename grid, unsigned int mom_p> void momenta_grid(grid & gd,typename grid::stype (& mom_tot)[grid::dims])
24{
25 auto it = gd.getDomainGhostIterator();
26
27 for (size_t i = 0 ; i < grid::dims ; i++)
28 mom_tot[i] = 0.0;
29
30 while (it.isNext())
31 {
32 auto key = it.get();
33 auto key_g = gd.getGKey(key);
34
35 for (size_t i = 0 ; i < grid::dims ; i++)
36 {
37 typename grid::stype coord = gd.spacing(i)*key_g.get(i);
38
39 mom_tot[i] += boost::math::pow<mom_p>(coord)*gd.template getProp<0>(key);
40 }
41
42 ++it;
43 }
44}
45
46template<typename grid, unsigned int mom_p> void momenta_grid_domain(grid & gd,typename grid::stype (& mom_tot)[grid::dims])
47{
48 auto it = gd.getDomainIterator();
49
50 for (size_t i = 0 ; i < grid::dims ; i++)
51 mom_tot[i] = 0.0;
52
53 while (it.isNext())
54 {
55 auto key = it.get();
56 auto key_g = gd.getGKey(key);
57
58 for (size_t i = 0 ; i < grid::dims ; i++)
59 {
60 typename grid::stype coord = gd.spacing(i)*key_g.get(i);
61
62 mom_tot[i] += boost::math::pow<mom_p>(coord)*gd.template getProp<0>(key);
63 }
64
65 ++it;
66 }
67}
68
69template<typename vector, unsigned int mom_p> void momenta_vector(vector & vd,typename vector::stype (& mom_tot)[vector::dims])
70{
71 auto it = vd.getDomainIterator();
72
73 for (size_t i = 0 ; i < vector::dims ; i++)
74 mom_tot[i] = 0.0;
75
76 while (it.isNext())
77 {
78 auto key = it.get();
79
80 for (size_t i = 0 ; i < vector::dims ; i++)
81 {
82 typename vector::stype coord = vd.getPos(key)[i];
83
84 mom_tot[i] += boost::math::pow<mom_p>(coord)*vd.template getProp<0>(key);
85 }
86
87 ++it;
88 }
89}
90
91template<unsigned int dim, typename T, typename grid, typename vector>
92void interp_test(grid & gd, vector & vd, bool single_particle)
93{
94 // Reset the grid
95
96 auto it2 = gd.getDomainGhostIterator();
97
98 while (it2.isNext())
99 {
100 auto key = it2.get();
101
102 gd.template get<0>(key) = 0.0;
103
104 ++it2;
105 }
106
107 interpolate<vector,grid,mp4_kernel<float>> inte(vd,gd);
108
109 if (single_particle == false)
110 {inte.template p2m<0,0>(vd,gd);}
111 else
112 {
113 auto it = vd.getDomainIterator();
114
115 while (it.isNext())
116 {
117 auto p = it.get();
118
119 inte.template p2m<0,0>(vd,gd,p);
120
121 ++it;
122 }
123 }
124
125 T mg[dim];
126 T mv[dim];
127
128 momenta_grid<grid,0>(gd,mg);
129 momenta_vector<vector,0>(vd,mv);
130
131 for (size_t i = 0 ; i < dim ; i++)
132 {BOOST_REQUIRE_CLOSE(mg[i],mv[i],0.001);}
133
134 momenta_grid<grid,1>(gd,mg);
135 momenta_vector<vector,1>(vd,mv);
136
137 for (size_t i = 0 ; i < dim ; i++)
138 {BOOST_REQUIRE_CLOSE(mg[i],mv[i],0.001);}
139
140 momenta_grid<grid,2>(gd,mg);
141 momenta_vector<vector,2>(vd,mv);
142
143 for (size_t i = 0 ; i < dim ; i++)
144 {BOOST_REQUIRE_CLOSE(mg[i],mv[i],0.001);}
145}
146
147BOOST_AUTO_TEST_CASE( interpolation_full_single_test_2D )
148{
149 Box<2,float> domain({0.0,0.0},{1.0,1.0});
150 size_t sz[2] = {64,64};
151
152 Ghost<2,long int> gg(2);
153 Ghost<2,float> gv(0.01);
154
155 size_t bc_v[2] = {PERIODIC,PERIODIC};
156
157 vector_dist<2,float,aggregate<float>> vd(65536,domain,bc_v,gv);
158 grid_dist_id<2,float,aggregate<float>> gd(vd.getDecomposition(),sz,gg);
159
160 // set one particle on vd
161
162 auto it = vd.getDomainIterator();
163
164 while (it.isNext())
165 {
166 auto p = it.get();
167
168 vd.getPos(p)[0] = (double)rand()/RAND_MAX;
169 vd.getPos(p)[1] = (double)rand()/RAND_MAX;
170
171 vd.getProp<0>(p) = 5.0;
172
173 ++it;
174 }
175
176 vd.map();
177
178 interp_test<2,float>(gd,vd,true);
179}
180
181
182BOOST_AUTO_TEST_CASE( interpolation_full_test_2D )
183{
184 Box<2,float> domain({0.0,0.0},{1.0,1.0});
185 size_t sz[2] = {64,64};
186
187 Ghost<2,long int> gg(2);
188 Ghost<2,float> gv(0.01);
189
190 size_t bc_v[2] = {PERIODIC,PERIODIC};
191
192 {
193 vector_dist<2,float,aggregate<float>> vd(4096,domain,bc_v,gv);
194 grid_dist_id<2,float,aggregate<float>> gd(vd.getDecomposition(),sz,gg);
195
196 // set one particle on vd
197
198 auto it = vd.getDomainIterator();
199
200 while (it.isNext())
201 {
202 auto p = it.get();
203
204 vd.getPos(p)[0] = (double)rand()/RAND_MAX;
205 vd.getPos(p)[1] = (double)rand()/RAND_MAX;
206
207 vd.getProp<0>(p) = 5.0;
208
209 ++it;
210 }
211
212 vd.map();
213
214 // Reset the grid
215 interp_test<2,float>(gd,vd,false);
216
217 float mg[2];
218 float mv[2];
219
220 auto & v_cl = create_vcluster();
221
222 interpolate<decltype(vd),decltype(gd),mp4_kernel<float>> inte(vd,gd);
223
224 // We have to do a ghost get before interpolating m2p
225 // Before doing mesh to particle particle must be arranged
226 // into a grid like
227
228 vd.clear();
229
230 auto it4 = vd.getGridIterator(sz);
231
232 while(it4.isNext())
233 {
234 auto key = it4.get();
235
236 vd.add();
237
238 vd.getLastPos()[0] = key.get(0) * it4.getSpacing(0) + domain.getLow(0) + 0.1*it4.getSpacing(0);
239 vd.getLastPos()[1] = key.get(1) * it4.getSpacing(1) + domain.getLow(1) + 0.1*it4.getSpacing(1);
240
241 vd.getLastProp<0>() = 0.0;
242
243 ++it4;
244 }
245
246 // Reset also the grid
247
248 auto it5 = gd.getDomainGhostIterator();
249
250 while(it5.isNext())
251 {
252 auto key = it5.get();
253
254 gd.get<0>(key) = 0.0;
255
256 ++it5;
257 }
258 gd.ghost_get<0>();
259
260 grid_key_dx<2> start({3,3});
261 grid_key_dx<2> stop({(long int)gd.size(0) - 4,(long int)gd.size(1) - 4});
262
263 auto it6 = gd.getSubDomainIterator(start,stop);
264 while(it6.isNext())
265 {
266 auto key = it6.get();
267
268 gd.get<0>(key) = 5.0/*(double)rand()/RAND_MAX;*/;
269
270 ++it6;
271 }
272 gd.ghost_get<0>();
273
274 vd.map();
275 gd.ghost_get<0>();
276 inte.m2p<0,0>(gd,vd);
277
278 momenta_grid_domain<decltype(gd),0>(gd,mg);
279 momenta_vector<decltype(vd),0>(vd,mv);
280
281 v_cl.sum(mg[0]);
282 v_cl.sum(mg[1]);
283 v_cl.sum(mv[0]);
284 v_cl.sum(mv[1]);
285 v_cl.execute();
286
287 BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
288 BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
289
290 momenta_grid_domain<decltype(gd),1>(gd,mg);
291 momenta_vector<decltype(vd),1>(vd,mv);
292
293 v_cl.sum(mg[0]);
294 v_cl.sum(mg[1]);
295 v_cl.sum(mv[0]);
296 v_cl.sum(mv[1]);
297 v_cl.execute();
298
299 BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
300 BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
301
302 momenta_grid_domain<decltype(gd),2>(gd,mg);
303 momenta_vector<decltype(vd),2>(vd,mv);
304
305 v_cl.sum(mg[0]);
306 v_cl.sum(mg[1]);
307 v_cl.sum(mv[0]);
308 v_cl.sum(mv[1]);
309 v_cl.execute();
310
311 BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
312 BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
313
314 }
315}
316
317BOOST_AUTO_TEST_CASE( interpolation_full_single_test_3D )
318{
319 Box<3,double> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
320 size_t sz[3] = {64,64,64};
321
322 Ghost<3,long int> gg(2);
323 Ghost<3,double> gv(0.01);
324
325 size_t bc_v[3] = {PERIODIC,PERIODIC,PERIODIC};
326
327 vector_dist<3,double,aggregate<double>> vd(65536,domain,bc_v,gv);
328 grid_dist_id<3,double,aggregate<double>> gd(vd.getDecomposition(),sz,gg);
329
330 // set one particle on vd
331
332 auto it = vd.getDomainIterator();
333
334 while (it.isNext())
335 {
336 auto p = it.get();
337
338 // coverty[dont_call]
339 vd.getPos(p)[0] = (double)rand()/RAND_MAX;
340 // coverty[dont_call]
341 vd.getPos(p)[1] = (double)rand()/RAND_MAX;
342 // coverty[dont_call]
343 vd.getPos(p)[2] = (double)rand()/RAND_MAX;
344
345 vd.getProp<0>(p) = 5.0;
346
347 ++it;
348 }
349
350 vd.map();
351
352 // Reset the grid
353 interp_test<3,double>(gd,vd,true);
354}
355
356BOOST_AUTO_TEST_CASE( interpolation_full_test_3D )
357{
358 Box<3,double> domain({0.0,0.0,0.0},{1.0,1.0,1.0});
359 size_t sz[3] = {64,64,64};
360
361 Ghost<3,long int> gg(2);
362 Ghost<3,double> gv(0.01);
363
364 size_t bc_v[3] = {PERIODIC,PERIODIC,PERIODIC};
365
366 {
367 vector_dist<3,double,aggregate<double>> vd(65536,domain,bc_v,gv);
368 grid_dist_id<3,double,aggregate<double>> gd(vd.getDecomposition(),sz,gg);
369
370 // set one particle on vd
371
372 auto it = vd.getDomainIterator();
373
374 while (it.isNext())
375 {
376 auto p = it.get();
377
378 vd.getPos(p)[0] = (double)rand()/RAND_MAX;
379 vd.getPos(p)[1] = (double)rand()/RAND_MAX;
380 vd.getPos(p)[2] = (double)rand()/RAND_MAX;
381
382 vd.getProp<0>(p) = 5.0;
383
384 ++it;
385 }
386
387 vd.map();
388
389 // Reset the grid
390
391 // Reset the grid
392 interp_test<3,double>(gd,vd,false);
393
394 auto & v_cl = create_vcluster();
395 double mg[3];
396 double mv[3];
397 interpolate<decltype(vd),decltype(gd),mp4_kernel<double>> inte(vd,gd);
398
399 // We have to do a ghost get before interpolating m2p
400 // Before doing mesh to particle particle must be arranged
401 // into a grid like
402
403 vd.clear();
404
405 auto it4 = vd.getGridIterator(sz);
406
407 while(it4.isNext())
408 {
409 auto key = it4.get();
410
411 vd.add();
412
413 vd.getLastPos()[0] = key.get(0) * it4.getSpacing(0) + domain.getLow(0) + 0.1*it4.getSpacing(0);
414 vd.getLastPos()[1] = key.get(1) * it4.getSpacing(1) + domain.getLow(1) + 0.1*it4.getSpacing(1);
415 vd.getLastPos()[2] = key.get(2) * it4.getSpacing(2) + domain.getLow(2) + 0.1*it4.getSpacing(2);
416
417 vd.getLastProp<0>() = 0.0;
418
419 ++it4;
420 }
421
422 // Reset also the grid
423
424 auto it5 = gd.getDomainGhostIterator();
425
426 while(it5.isNext())
427 {
428 auto key = it5.get();
429
430 gd.get<0>(key) = 0.0;
431
432 ++it5;
433 }
434 gd.ghost_get<0>();
435
436 grid_key_dx<3> start({3,3,3});
437 grid_key_dx<3> stop({(long int)gd.size(0) - 4,(long int)gd.size(1) - 4,(long int)gd.size(2) - 4});
438
439 auto it6 = gd.getSubDomainIterator(start,stop);
440 while(it6.isNext())
441 {
442 auto key = it6.get();
443
444 gd.get<0>(key) = 5.0;
445
446 ++it6;
447 }
448 gd.ghost_get<0>();
449
450 vd.map();
451 gd.ghost_get<0>();
452 inte.m2p<0,0>(gd,vd);
453
454 momenta_grid_domain<decltype(gd),0>(gd,mg);
455 momenta_vector<decltype(vd),0>(vd,mv);
456
457 v_cl.sum(mg[0]);
458 v_cl.sum(mg[1]);
459 v_cl.sum(mg[2]);
460 v_cl.sum(mv[0]);
461 v_cl.sum(mv[1]);
462 v_cl.sum(mv[2]);
463 v_cl.execute();
464
465 BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
466 BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
467 BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
468
469 momenta_grid_domain<decltype(gd),1>(gd,mg);
470 momenta_vector<decltype(vd),1>(vd,mv);
471
472 v_cl.sum(mg[0]);
473 v_cl.sum(mg[1]);
474 v_cl.sum(mg[2]);
475 v_cl.sum(mv[0]);
476 v_cl.sum(mv[1]);
477 v_cl.sum(mv[2]);
478 v_cl.execute();
479
480 BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
481 BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
482 BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
483
484 momenta_grid_domain<decltype(gd),2>(gd,mg);
485 momenta_vector<decltype(vd),2>(vd,mv);
486
487 v_cl.sum(mg[0]);
488 v_cl.sum(mg[1]);
489 v_cl.sum(mg[2]);
490 v_cl.sum(mv[0]);
491 v_cl.sum(mv[1]);
492 v_cl.sum(mv[2]);
493 v_cl.execute();
494
495 BOOST_REQUIRE_CLOSE(mg[0],mv[0],0.001);
496 BOOST_REQUIRE_CLOSE(mg[1],mv[1],0.001);
497 BOOST_REQUIRE_CLOSE(mg[2],mv[2],0.001);
498
499 }
500}
501
502BOOST_AUTO_TEST_CASE( int_kernel_test )
503{
504 mp4_kernel<float> mp4;
505
506 float tot = 0.0;
507
508 // Check momenta 0
509
510 tot += mp4.value(-1.3f,0);
511 tot += mp4.value(-0.3f,1);
512 tot += mp4.value(0.7f,2);
513 tot += mp4.value(1.7f,3);
514
515 BOOST_REQUIRE_CLOSE(tot,1.0f,0.001);
516
517 // Check momenta 1
518
519 tot = 0.0;
520
521 tot += -1.3f*mp4.value(-1.3f,0);
522 tot += -0.3f*mp4.value(-0.3f,1);
523 tot += 0.7f*mp4.value(0.7f,2);
524 tot += 1.7f*mp4.value(1.7f,3);
525
526 BOOST_REQUIRE_SMALL(tot,0.001f);
527
528 // Check momenta 2
529
530 tot = 0.0;
531
532 tot += (1.3f)*(1.3f)*mp4.value(-1.3f,0);
533 tot += (0.3f)*(0.3f)*mp4.value(-0.3f,1);
534 tot += (0.7f)*(0.7f)*mp4.value(0.7f,2);
535 tot += (1.7f)*(1.7f)*mp4.value(1.7f,3);
536
537 BOOST_REQUIRE_SMALL(tot,0.001f);
538
539
540 //////// Check zeta 1
541
542 tot = 0.0;
543
544 z_kernel<float,1> zk1;
545
546 tot += zk1.value(-0.3f,0);
547 tot += zk1.value(0.7f,1);
548
549 BOOST_REQUIRE_CLOSE(tot,1.0f,0.001);
550
551 //////// zeta 2 is equivalent to mp4 we do not test
552
553 //////// zeta 3
554
555 z_kernel<float,3> zk3;
556
557 tot = 0.0;
558
559 // Check momenta 0
560
561 tot += zk3.value(-2.3f,0);
562 tot += zk3.value(-1.3f,1);
563 tot += zk3.value(-0.3f,2);
564 tot += zk3.value(0.7f,3);
565 tot += zk3.value(1.7f,4);
566 tot += zk3.value(2.7f,5);
567
568 BOOST_REQUIRE_CLOSE(tot,1.0f,0.001);
569
570 // Check momenta 1
571
572 tot = 0.0;
573
574 tot += -2.3*zk3.value(-2.3f,0);
575 tot += -1.3*zk3.value(-1.3f,1);
576 tot += -0.3*zk3.value(-0.3f,2);
577 tot += 0.7*zk3.value(0.7f,3);
578 tot += 1.7*zk3.value(1.7f,4);
579 tot += 2.7*zk3.value(2.7f,5);
580
581 BOOST_REQUIRE_SMALL(tot,0.001f);
582
583 // Check momenta 2
584
585 tot = 0.0;
586
587 tot += 2.3*2.3*zk3.value(-2.3f,0);
588 tot += 1.3*1.3*zk3.value(-1.3f,1);
589 tot += 0.3*0.3*zk3.value(-0.3f,2);
590 tot += 0.7*0.7*zk3.value(0.7f,3);
591 tot += 1.7*1.7*zk3.value(1.7f,4);
592 tot += 2.7*2.7*zk3.value(2.7f,5);
593
594 BOOST_REQUIRE_SMALL(tot,0.001f);
595
596 // Check momenta 3
597
598 tot = 0.0;
599
600 tot += -2.3*-2.3*-2.3*zk3.value(-2.3f,0);
601 tot += -1.3*-1.3*-1.3*zk3.value(-1.3f,1);
602 tot += -0.3*-0.3*-0.3*zk3.value(-0.3f,2);
603 tot += 0.7*0.7*0.7*zk3.value(0.7f,3);
604 tot += 1.7*1.7*1.7*zk3.value(1.7f,4);
605 tot += 2.7*2.7*2.7*zk3.value(2.7f,5);
606
607 BOOST_REQUIRE_SMALL(tot,0.001f);
608
609
610 // z4
611
612 z_kernel<float,4> zk4;
613
614 // Check momenta 0
615
616 tot = 0.0;
617
618 tot += zk4.value(-3.3f,0);
619 tot += zk4.value(-2.3f,1);
620 tot += zk4.value(-1.3f,2);
621 tot += zk4.value(-0.3f,3);
622 tot += zk4.value(0.7f,4);
623 tot += zk4.value(1.7f,5);
624 tot += zk4.value(2.7f,6);
625 tot += zk4.value(3.7f,7);
626
627 BOOST_REQUIRE_CLOSE(tot,1.0f,0.001);
628
629 // Check momenta 1
630
631 tot = 0.0;
632
633 tot += -3.3*zk4.value(-3.3f,0);
634 tot += -2.3*zk4.value(-2.3f,1);
635 tot += -1.3*zk4.value(-1.3f,2);
636 tot += -0.3*zk4.value(-0.3f,3);
637 tot += 0.7*zk4.value(0.7f,4);
638 tot += 1.7*zk4.value(1.7f,5);
639 tot += 2.7*zk4.value(2.7f,6);
640 tot += 3.7*zk4.value(3.7f,7);
641
642 BOOST_REQUIRE_SMALL(tot,0.001f);
643
644 // Check momenta 2
645
646 tot = 0.0;
647
648 tot += 3.3*3.3*zk4.value(-3.3f,0);
649 tot += 2.3*2.3*zk4.value(-2.3f,1);
650 tot += 1.3*1.3*zk4.value(-1.3f,2);
651 tot += 0.3*0.3*zk4.value(-0.3f,3);
652 tot += 0.7*0.7*zk4.value(0.7f,4);
653 tot += 1.7*1.7*zk4.value(1.7f,5);
654 tot += 2.7*2.7*zk4.value(2.7f,6);
655 tot += 3.7*3.7*zk4.value(3.7f,7);
656
657 BOOST_REQUIRE_SMALL(tot,0.001f);
658
659 // Check momenta 3
660
661 tot = 0.0;
662
663 tot += -3.3*-3.3*-3.3*zk4.value(-3.3f,0);
664 tot += -2.3*-2.3*-2.3*zk4.value(-2.3f,1);
665 tot += -1.3*-1.3*-1.3*zk4.value(-1.3f,2);
666 tot += -0.3*-0.3*-0.3*zk4.value(-0.3f,3);
667 tot += 0.7*0.7*0.7*zk4.value(0.7f,4);
668 tot += 1.7*1.7*1.7*zk4.value(1.7f,5);
669 tot += 2.7*2.7*2.7*zk4.value(2.7f,6);
670 tot += 3.7*3.7*3.7*zk4.value(3.7f,7);
671
672 BOOST_REQUIRE_SMALL(tot,0.001f);
673
674 // Check momenta 4
675
676 tot = 0.0;
677
678 tot += -3.3*-3.3*-3.3*-3.3*zk4.value(-3.3f,0);
679 tot += -2.3*-2.3*-2.3*-2.3*zk4.value(-2.3f,1);
680 tot += -1.3*-1.3*-1.3*-1.3*zk4.value(-1.3f,2);
681 tot += -0.3*-0.3*-0.3*-0.3*zk4.value(-0.3f,3);
682 tot += 0.7*0.7*0.7*0.7*zk4.value(0.7f,4);
683 tot += 1.7*1.7*1.7*1.7*zk4.value(1.7f,5);
684 tot += 2.7*2.7*2.7*2.7*zk4.value(2.7f,6);
685 tot += 3.7*3.7*3.7*3.7*zk4.value(3.7f,7);
686
687 BOOST_REQUIRE_SMALL(tot,0.001f);
688}
689
690BOOST_AUTO_TEST_SUITE_END()
691
692
693#endif /* OPENFPM_NUMERICS_SRC_INTERPOLATION_INTERPOLATION_UNIT_TESTS_HPP_ */
694