OpenFPM  5.2.0
Project that contain the implementation of distributed structures
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 
46 struct no_equation
47 {};
48 
49 BOOST_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 
516 BOOST_AUTO_TEST_SUITE_END()
Box
This class represent an N-dimensional box.
Definition: Box.hpp:60
Ghost
Definition: Ghost.hpp:40
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:279
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:204
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