OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
FD_op_Tests.cpp
1//
2// Created by Abhinav Singh on 01.05.20.
3//
4
5
6#include "config.h"
7
8#define BOOST_TEST_DYN_LINK
9
10#include "util/util_debug.hpp"
11#include <boost/test/unit_test.hpp>
12#include <iostream>
13#include "Grid/grid_dist_id.hpp"
14#include "FD_op.hpp"
15#include "Grid/staggered_dist_grid.hpp"
16//#include "FDScheme.hpp"
17//#include "Solvers/petsc_solver.hpp"
18//#include "DCPSE_op/EqnsStruct.hpp"
19
21/*struct equations2d1 {
22
24 static const unsigned int dims=2;
26 static const unsigned int nvar=1;
27
29 static const bool boundary[];
30
32 typedef double stype;
33
35 typedef grid_dist_id<2, double, aggregate<double, double, double>> b_part;
36
38 typedef SparseMatrix<double, int, PETSC_BASE> SparseMatrix_type;
39
41 typedef Vector<double, PETSC_BASE> Vector_type;
42
43 typedef petsc_solver<double> solver_type;
44};*/
45
46struct no_equation
47{};
48
49BOOST_AUTO_TEST_SUITE(fd_op_suite_tests)
50
51 BOOST_AUTO_TEST_CASE(fd_op_tests) {
52 size_t edgeSemiSize = 80;
53 const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
54 Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
55 periodicity<2> bc({NON_PERIODIC, NON_PERIODIC});
56 double spacing[2];
57 spacing[0] = 2 * M_PI / (sz[0] - 1);
58 spacing[1] = 2 * M_PI / (sz[1] - 1);
59 Ghost<2, long int> ghost(1);
60
61 //std::cout << "Spacing: " << spacing[0] << " " << spacing[1] << std::endl;
62
63 grid_dist_id<2, double, aggregate<double, double, double>> domain(sz, box,ghost,bc);
64
65 BOOST_TEST_MESSAGE("Init domain...");
66 auto it = domain.getDomainIterator();
67 while (it.isNext())
68 {
69 auto key_l = it.get();
70 auto key = it.getGKey(key_l);
71 mem_id i = key.get(0);
72 double x = i * spacing[0];
73 mem_id j = key.get(1);
74 double y = j * spacing[1];
75 // Here fill the function value P
76 domain.template getProp<0>(key_l) = sin(x) + sin(y);
77 domain.template getProp<1>(key_l) = 0;
78 // Here fill the validation value for Df/Dx in property 3
79
80 domain.template getProp<2>(key_l) = cos(x) + cos(y) + sin(x) + sin(y) + 5.0;
81
82 ++it;
83 }
84
85 domain.ghost_get<0>();
86
87 FD::Derivative_x Dx;
88 FD::Derivative_y Dy;
89
90 auto v = FD::getV<1>(domain);
91 auto P = FD::getV<0>(domain);
92
93 v = Dx(P) + Dy(P) + P + 5;
94
95 auto it2 = domain.getDomainIterator();
96
97 double worst = 0.0;
98
99 while (it2.isNext()) {
100 auto p = it2.get();
101
102 if (fabs(domain.getProp<1>(p) - domain.getProp<2>(p)) > worst) {
103 worst = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
104 }
105
106 ++it2;
107 }
108
109 //std::cout << "Maximum Error: " << worst << std::endl;
110
111 domain.write("g");
112
113 BOOST_REQUIRE(worst < 0.003);
114 }
115
116
117 BOOST_AUTO_TEST_CASE(fd_op_tests_vec_mat) {
118 size_t edgeSemiSize = 80;
119 const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
120 Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
121 periodicity<2> bc({NON_PERIODIC, NON_PERIODIC});
122 double spacing[2];
123 spacing[0] = 2 * M_PI / (sz[0] - 1);
124 spacing[1] = 2 * M_PI / (sz[1] - 1);
125 Ghost<2, long int> ghost(1);
126
127 //std::cout << "Spacing: " << spacing[0] << " " << spacing[1] << std::endl;
128
129 grid_dist_id<2, double, aggregate<double, double, double,double[2],double[2][2],double[2][2][2]>> domain(sz, box,ghost,bc);
130
131 BOOST_TEST_MESSAGE("Init domain...");
132 auto it = domain.getDomainIterator();
133 while (it.isNext())
134 {
135 auto key_l = it.get();
136 auto key = it.getGKey(key_l);
137 mem_id i = key.get(0);
138 double x = i * spacing[0];
139 mem_id j = key.get(1);
140 double y = j * spacing[1];
141 // Here fill the function value P
142 domain.template getProp<3>(key_l)[0] = sin(x);
143 domain.template getProp<3>(key_l)[1] = sin(x);
144 domain.template getProp<1>(key_l) = 0;
145 // Here fill the validation value for Df/Dx in property 3
146
147 domain.template getProp<2>(key_l) = cos(x);
148
149 ++it;
150 }
151
152 domain.ghost_get<0,3>();
153
154 FD::Derivative_x Dx;
155 FD::Derivative_y Dy;
156
157 auto v = FD::getV<1>(domain);
158 auto P = FD::getV<0>(domain);
159
160 auto vec = FD::getV<3>(domain);
161 auto Mat = FD::getV<4>(domain);
162 auto Mat3 = FD::getV<5>(domain);
163
164 Mat[0][1] = Dx(vec[0]);
165
166 Mat[1][0] = vec[0];
167
168 domain.ghost_get<4>();
169
170 Mat[0][0] = Dx(Mat[1][0]);
171 Mat3[0][0][0] = Dx(vec[0]);
172 Mat3[0][1][0] = Dx(Mat[1][0]);
173
174 auto it2 = domain.getDomainIterator();
175
176 double worst = 0.0;
177
178 while (it2.isNext())
179 {
180 auto p = it2.get();
181
182 if (fabs(domain.getProp<4>(p)[0][1] - domain.getProp<2>(p)) > worst)
183 {
184 worst = fabs(domain.getProp<4>(p)[0][1] - domain.getProp<2>(p));
185 }
186
187 if (fabs(domain.getProp<4>(p)[1][0] - domain.getProp<3>(p)[0]) > worst)
188 {
189 worst = fabs(domain.getProp<4>(p)[1][0] - domain.getProp<3>(p)[0]);
190 }
191
192 if (fabs(domain.getProp<4>(p)[0][0] - domain.getProp<2>(p)) > worst)
193 {
194 worst = fabs(domain.getProp<4>(p)[0][0] - domain.getProp<2>(p));
195 }
196
197
199
200 if (fabs(domain.getProp<5>(p)[0][0][0] - domain.getProp<2>(p)) > worst)
201 {
202 worst = fabs(domain.getProp<5>(p)[0][0][0] - domain.getProp<2>(p));
203 }
204
205 if (fabs(domain.getProp<5>(p)[0][1][0] - domain.getProp<2>(p)) > worst)
206 {
207 worst = fabs(domain.getProp<5>(p)[0][1][0] - domain.getProp<2>(p));
208 }
209
210 ++it2;
211 }
212
213 BOOST_REQUIRE(worst < 0.003);
214 }
215
216 BOOST_AUTO_TEST_CASE(lalpacian_test) {
217 size_t edgeSemiSize = 80;
218 const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
219 Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
220 periodicity<2> bc({NON_PERIODIC, NON_PERIODIC});
221 double spacing[2];
222 spacing[0] = 2 * M_PI / (sz[0] - 1);
223 spacing[1] = 2 * M_PI / (sz[1] - 1);
224 Ghost<2, long int> ghost(1);
225
226 //std::cout << "Spacing: " << spacing[0] << " " << spacing[1] << std::endl;
227
228 grid_dist_id<2, double, aggregate<double, double, double>> domain(sz, box,ghost,bc);
229
230 BOOST_TEST_MESSAGE("Init domain...");
231 auto it = domain.getDomainIterator();
232 while (it.isNext())
233 {
234 auto key_l = it.get();
235 auto key = it.getGKey(key_l);
236 mem_id i = key.get(0);
237 double x = M_PI / 2 + i * spacing[0];
238 mem_id j = M_PI / 2 + key.get(1);
239 double y = j * spacing[1];
240 // Here fill the function value P
241 domain.template getProp<0>(key_l) = sin(x) + sin(y);
242 domain.template getProp<1>(key_l) = 0;
243 // Here fill the validation value for Df/Dx in property 3
244
245 domain.template getProp<2>(key_l) = -sin(x) -sin(y) + cos(x) + cos(y) + 5.0;
246
247 ++it;
248 }
249
250 domain.ghost_get<0>();
251
252 FD::Derivative_x Dx;
253 FD::Derivative_y Dy;
254 FD::Lap L;
255
256 auto v = FD::getV<1>(domain);
257 auto P = FD::getV<0>(domain);
258
259 v = L(P) + Dx(P) + Dy(P) + 5;
260
261 auto it2 = domain.getDomainIterator();
262
263 double worst = 0.0;
264
265 while (it2.isNext()) {
266 auto p = it2.get();
267
268 if (fabs(domain.getProp<1>(p) - domain.getProp<2>(p)) > worst) {
269 worst = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
270 }
271
272 ++it2;
273 }
274
275 //std::cout << "Maximum Error: " << worst << std::endl;
276
277 //domain.write("g");
278
279 BOOST_REQUIRE(worst < 0.003);
280 }
281
282
283 BOOST_AUTO_TEST_CASE(Slice_test) {
284 size_t edgeSemiSize = 80;
285 const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
286 Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
287 periodicity<2> bc({NON_PERIODIC, NON_PERIODIC});
288 double spacing[2];
289 spacing[0] = 2 * M_PI / (sz[0] - 1);
290 spacing[1] = 2 * M_PI / (sz[1] - 1);
291 Ghost<2, long int> ghost(1);
292
293 //std::cout << "Spacing: " << spacing[0] << " " << spacing[1] << std::endl;
294
295 grid_dist_id<2, double, aggregate<double[2], double[2], double[2],double>> domain(sz, box,ghost,bc);
296
297 BOOST_TEST_MESSAGE("Init domain...");
298 auto it = domain.getDomainIterator();
299 while (it.isNext())
300 {
301 auto key_l = it.get();
302 auto key = it.getGKey(key_l);
303 mem_id i = key.get(0);
304 double x = M_PI / 2 + i * spacing[0];
305 mem_id j = M_PI / 2 + key.get(1);
306 double y = j * spacing[1];
307 // Here fill the function value P
308 domain.template getProp<0>(key_l)[0] = sin(x) + sin(y);
309 domain.template getProp<0>(key_l)[1] = cos(x) + cos(y);
310 domain.template getProp<2>(key_l)[0] = 0;
311 domain.template getProp<2>(key_l)[1] = 0;
312 domain.template getProp<3>(key_l) = 0;
313 // Here fill the validation value for Df/Dx in property 3
314
315 domain.template getProp<1>(key_l)[0] = -sin(x) -sin(y) + cos(x) + cos(y) + 5.0;
316 domain.template getProp<1>(key_l)[1] = -cos(x) -cos(y) - sin(x) - sin(y) + 5.0;
317
318 ++it;
319 }
320
321 domain.ghost_get<0>();
322
323 FD::Derivative_x Dx;
324 FD::Derivative_y Dy;
325 FD::Lap L;
326
327
328 auto V1 = FD::getV<0>(domain);
329 auto V2 = FD::getV<2>(domain);
330 auto P = FD::getV<3>(domain);
331 constexpr int x=0,y=1;
332 //V2[x] = Dx(V1[y]) + Dy(V1[y]) + 5;
333 V2[x] = L(V1[y]) + Dx(V1[y]) + Dy(V1[y]) + 5;
334 V2[y] = L(V1[x]) + Dx(V1[x]) + Dy(V1[x]) + 5;
335
336 auto it2 = domain.getDomainIterator();
337
338 double worst = 0.0;
339
340 while (it2.isNext()) {
341 auto p = it2.get();
342
343 if (fabs(domain.getProp<1>(p)[0] - domain.getProp<2>(p)[1]) > worst) {
344 worst = fabs(domain.getProp<1>(p)[0] - domain.getProp<2>(p)[1]);
345 }
346
347 ++it2;
348 }
349
350 //std::cout << "Maximum Error: " << worst << std::endl;
351
352 //domain.write("g");
353
354 BOOST_REQUIRE(worst < 0.003);
355 BOOST_REQUIRE(worst != 0);
356
357 }
358
359 BOOST_AUTO_TEST_CASE(fd_op_tests_staggered) {
360
361 constexpr int phi_ = 0;
362 constexpr int ux_ = 1;
363
364 size_t edgeSemiSize = 80;
365 const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
366 Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
367 periodicity<2> bc({NON_PERIODIC, NON_PERIODIC});
368 double spacing[2];
369 spacing[0] = 2 * M_PI / (sz[0] - 1);
370 spacing[1] = 2 * M_PI / (sz[1] - 1);
371 Ghost<2, long int> ghost(1);
372
373 //std::cout << "Spacing: " << spacing[0] << " " << spacing[1] << std::endl;
374
375 staggered_grid_dist<2, double, aggregate<double, double>> domain(sz, box,ghost,bc);
376
377 comb<2> c({(char)-1,(char)-1});
378 openfpm::vector<comb<2>> lc;
379 lc.add(c);
380
381 domain.setStagPosition<phi_>(lc);
382
383 comb<2> c2({(char)0,(char)-1});
384 lc.clear();
385 lc.add(c2);
386
387 domain.setStagPosition<ux_>(lc);
388
389 BOOST_TEST_MESSAGE("Init domain...");
390 auto it = domain.getDomainIterator();
391 while (it.isNext())
392 {
393 auto key_l = it.get();
394 auto key = it.getGKey(key_l);
395 mem_id i = key.get(0);
396 mem_id j = key.get(1);
397 // Here fill the function value P
398 domain.template getProp<phi_>(key_l) = j % 2;
399
400 ++it;
401 }
402
403 domain.ghost_get<0>();
404
405 FD::Derivative_x Dx;
406 FD::Derivative_y Dy;
407
408 auto ux = FD::getV_stag<ux_>(domain);
409 auto phi = FD::getV_stag<phi_>(domain);
410
411 grid_key_dx<2> start({0ul,0u});
412 grid_key_dx<2> stop({(long int)(sz[0]-2),(long int)(sz[1]-2)});
413
414 auto it2 = domain.getSubDomainIterator(start,stop);
415 while (it2.isNext())
416 {
417 auto key_l = it2.get();
418
419 double test = phi.value(key_l,c2);
420
421 BOOST_REQUIRE_EQUAL(test,0.5);
422
423 ++it2;
424 }
425 }
426
427 BOOST_AUTO_TEST_CASE(fd_op_tests_derivative_staggered) {
428
429 constexpr int x = 0;
430 constexpr int y = 1;
431
432 constexpr int phi_ = 0;
433 constexpr int ux_ = 1;
434
435 size_t edgeSemiSize = 80;
436 const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
437 Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
438 periodicity<2> bc({NON_PERIODIC, NON_PERIODIC});
439 double spacing[2];
440 spacing[0] = 2 * M_PI / (sz[0] - 1);
441 spacing[1] = 2 * M_PI / (sz[1] - 1);
442 Ghost<2, long int> ghost(1);
443
444 //std::cout << "Spacing: " << spacing[0] << " " << spacing[1] << std::endl;
445
446 staggered_grid_dist<2, double, aggregate<double, double>> domain(sz, box,ghost,bc);
447
448 comb<2> c({(char)-1,(char)-1});
449 openfpm::vector<comb<2>> lc;
450 lc.add(c);
451
452 domain.setStagPosition<phi_>(lc);
453
454 comb<2> c2({(char)0,(char)-1});
455 lc.clear();
456 lc.add(c2);
457
458 domain.setStagPosition<ux_>(lc);
459
460 BOOST_TEST_MESSAGE("Init domain...");
461 auto it = domain.getDomainIterator();
462 while (it.isNext())
463 {
464 auto key_l = it.get();
465 auto key = it.getGKey(key_l);
466 mem_id i = key.get(0);
467 double x = i * spacing[0];
468 mem_id j = key.get(1);
469 double y = j * spacing[1];
470 // Here fill the function value P
471 domain.template getProp<phi_>(key_l) = x*x + 3.0*y*y;
472
473 ++it;
474 }
475
476 domain.ghost_get<0>();
477
478 FD::Derivative_x Dx;
479 FD::Derivative_y Dy;
480
481 auto ux = FD::getV_stag<ux_>(domain);
482 auto phi = FD::getV_stag<phi_>(domain);
483
484 ux = Dx(phi) + Dy(phi) + phi + 5;
485
486 grid_key_dx<2> start({1ul,1u});
487 grid_key_dx<2> stop({(long int)(sz[0]-2),(long int)(sz[1]-2)});
488
489 double worst_error = 0.0;
490
491 auto it2 = domain.getSubDomainIterator(start,stop);
492 while (it2.isNext())
493 {
494 auto key = it2.get();
495
496 double test = domain.template getProp<ux_>(key);
497
498 double dx = 0.5* (domain.template getProp<phi_>(key.move(x,1).move(y,1)) + domain.template getProp<phi_>(key.move(x,1)) - (domain.template getProp<phi_>(key.move(x,-1).move(y,1)) + domain.template getProp<phi_>(key.move(x,-1))) );
499 double dy = 0.5* (domain.template getProp<phi_>(key.move(y,1).move(y,1)) + domain.template getProp<phi_>(key.move(y,1)) - (domain.template getProp<phi_>(key.move(y,-1).move(y,1)) + domain.template getProp<phi_>(key.move(y,-1))) );
500
501 dx *= 0.5/spacing[0];
502 dy *= 0.5/spacing[1];
503
504 double phi_inte = 0.5*(domain.template getProp<phi_>(key.move(y,1)) + domain.template getProp<phi_>(key));
505
506 if (fabs(test - (dx + dy + phi_inte + 5.0)) > worst_error )
507 {worst_error = fabs(test - (dx + dy + phi_inte + 5.0));}
508
509 ++it2;
510 }
511
512 BOOST_REQUIRE(worst_error < 5e-13);
513 }
514
515
516BOOST_AUTO_TEST_SUITE_END()
Box
This class represent an N-dimensional box.
Definition Box.hpp:61
Point_test
Test structure used for several test.
Definition Point_test.hpp:106
grid_dist_id
This is a distributed grid.
Definition grid_dist_id.hpp:278
grid_key_dx
grid_key_dx is the key to access any element in the grid
Definition grid_key.hpp:19
grid_key_dx::get
__device__ __host__ index_type get(index_type i) const
Get the i index.
Definition grid_key.hpp:503
openfpm::vector
Implementation of 1-D std::vector like structure.
Definition map_vector.hpp:203
staggered_grid_dist
Implementation of the staggered grid.
Definition staggered_dist_grid.hpp:64
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
comb
Position of the element of dimension d in the hyper-cube of dimension dim.
Definition comb.hpp:35
no_equation
Specify the general characteristic of system to solve.
Definition FD_op_Tests.cpp:47
periodicity
Boundary conditions.
Definition common.hpp:22