OpenFPM_pdata  4.1.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(lalpacian_test) {
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>> 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 = M_PI / 2 + i * spacing[0];
139  mem_id j = M_PI / 2 + key.get(1);
140  double y = j * spacing[1];
141  // Here fill the function value P
142  domain.template getProp<0>(key_l) = sin(x) + sin(y);
143  domain.template getProp<1>(key_l) = 0;
144  // Here fill the validation value for Df/Dx in property 3
145 
146  domain.template getProp<2>(key_l) = -sin(x) -sin(y) + cos(x) + cos(y) + 5.0;
147 
148  ++it;
149  }
150 
151  domain.ghost_get<0>();
152 
153  FD::Derivative_x Dx;
154  FD::Derivative_y Dy;
155  FD::Lap L;
156 
157  auto v = FD::getV<1>(domain);
158  auto P = FD::getV<0>(domain);
159 
160  v = L(P) + Dx(P) + Dy(P) + 5;
161 
162  auto it2 = domain.getDomainIterator();
163 
164  double worst = 0.0;
165 
166  while (it2.isNext()) {
167  auto p = it2.get();
168 
169  if (fabs(domain.getProp<1>(p) - domain.getProp<2>(p)) > worst) {
170  worst = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
171  }
172 
173  ++it2;
174  }
175 
176  //std::cout << "Maximum Error: " << worst << std::endl;
177 
178  //domain.write("g");
179 
180  BOOST_REQUIRE(worst < 0.003);
181  }
182 
183 
184  BOOST_AUTO_TEST_CASE(Slice_test) {
185  size_t edgeSemiSize = 80;
186  const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
187  Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
188  periodicity<2> bc({NON_PERIODIC, NON_PERIODIC});
189  double spacing[2];
190  spacing[0] = 2 * M_PI / (sz[0] - 1);
191  spacing[1] = 2 * M_PI / (sz[1] - 1);
192  Ghost<2, long int> ghost(1);
193 
194  //std::cout << "Spacing: " << spacing[0] << " " << spacing[1] << std::endl;
195 
196  grid_dist_id<2, double, aggregate<double[2], double[2], double[2],double>> domain(sz, box,ghost,bc);
197 
198  BOOST_TEST_MESSAGE("Init domain...");
199  auto it = domain.getDomainIterator();
200  while (it.isNext())
201  {
202  auto key_l = it.get();
203  auto key = it.getGKey(key_l);
204  mem_id i = key.get(0);
205  double x = M_PI / 2 + i * spacing[0];
206  mem_id j = M_PI / 2 + key.get(1);
207  double y = j * spacing[1];
208  // Here fill the function value P
209  domain.template getProp<0>(key_l)[0] = sin(x) + sin(y);
210  domain.template getProp<0>(key_l)[1] = cos(x) + cos(y);
211  domain.template getProp<2>(key_l)[0] = 0;
212  domain.template getProp<2>(key_l)[1] = 0;
213  domain.template getProp<3>(key_l) = 0;
214  // Here fill the validation value for Df/Dx in property 3
215 
216  domain.template getProp<1>(key_l)[0] = -sin(x) -sin(y) + cos(x) + cos(y) + 5.0;
217  domain.template getProp<1>(key_l)[1] = -cos(x) -cos(y) - sin(x) - sin(y) + 5.0;
218 
219  ++it;
220  }
221 
222  domain.ghost_get<0>();
223 
224  FD::Derivative_x Dx;
225  FD::Derivative_y Dy;
226  FD::Lap L;
227 
228 
229  auto V1 = FD::getV<0>(domain);
230  auto V2 = FD::getV<2>(domain);
231  auto P = FD::getV<3>(domain);
232  constexpr int x=0,y=1;
233  //V2[x] = Dx(V1[y]) + Dy(V1[y]) + 5;
234  V2[x] = L(V1[y]) + Dx(V1[y]) + Dy(V1[y]) + 5;
235  V2[y] = L(V1[x]) + Dx(V1[x]) + Dy(V1[x]) + 5;
236 
237  auto it2 = domain.getDomainIterator();
238 
239  double worst = 0.0;
240 
241  while (it2.isNext()) {
242  auto p = it2.get();
243 
244  if (fabs(domain.getProp<1>(p)[0] - domain.getProp<2>(p)[1]) > worst) {
245  worst = fabs(domain.getProp<1>(p)[0] - domain.getProp<2>(p)[1]);
246  }
247 
248  ++it2;
249  }
250 
251  //std::cout << "Maximum Error: " << worst << std::endl;
252 
253  //domain.write("g");
254 
255  BOOST_REQUIRE(worst < 0.003);
256  BOOST_REQUIRE(worst != 0);
257 
258  }
259 
260  BOOST_AUTO_TEST_CASE(fd_op_tests_staggered) {
261 
262  constexpr int phi_ = 0;
263  constexpr int ux_ = 1;
264 
265  size_t edgeSemiSize = 80;
266  const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
267  Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
268  periodicity<2> bc({NON_PERIODIC, NON_PERIODIC});
269  double spacing[2];
270  spacing[0] = 2 * M_PI / (sz[0] - 1);
271  spacing[1] = 2 * M_PI / (sz[1] - 1);
272  Ghost<2, long int> ghost(1);
273 
274  //std::cout << "Spacing: " << spacing[0] << " " << spacing[1] << std::endl;
275 
276  staggered_grid_dist<2, double, aggregate<double, double>> domain(sz, box,ghost,bc);
277 
278  comb<2> c({(char)-1,(char)-1});
279  openfpm::vector<comb<2>> lc;
280  lc.add(c);
281 
282  domain.setStagPosition<phi_>(lc);
283 
284  comb<2> c2({(char)0,(char)-1});
285  lc.clear();
286  lc.add(c2);
287 
288  domain.setStagPosition<ux_>(lc);
289 
290  BOOST_TEST_MESSAGE("Init domain...");
291  auto it = domain.getDomainIterator();
292  while (it.isNext())
293  {
294  auto key_l = it.get();
295  auto key = it.getGKey(key_l);
296  mem_id i = key.get(0);
297  mem_id j = key.get(1);
298  // Here fill the function value P
299  domain.template getProp<phi_>(key_l) = j % 2;
300 
301  ++it;
302  }
303 
304  domain.ghost_get<0>();
305 
306  FD::Derivative_x Dx;
307  FD::Derivative_y Dy;
308 
309  auto ux = FD::getV_stag<ux_>(domain);
310  auto phi = FD::getV_stag<phi_>(domain);
311 
312  grid_key_dx<2> start({0ul,0u});
313  grid_key_dx<2> stop({(long int)(sz[0]-2),(long int)(sz[1]-2)});
314 
315  auto it2 = domain.getSubDomainIterator(start,stop);
316  while (it2.isNext())
317  {
318  auto key_l = it2.get();
319 
320  double test = phi.value(key_l,c2);
321 
322  BOOST_REQUIRE_EQUAL(test,0.5);
323 
324  ++it2;
325  }
326  }
327 
328  BOOST_AUTO_TEST_CASE(fd_op_tests_derivative_staggered) {
329 
330  constexpr int x = 0;
331  constexpr int y = 1;
332 
333  constexpr int phi_ = 0;
334  constexpr int ux_ = 1;
335 
336  size_t edgeSemiSize = 80;
337  const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
338  Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
339  periodicity<2> bc({NON_PERIODIC, NON_PERIODIC});
340  double spacing[2];
341  spacing[0] = 2 * M_PI / (sz[0] - 1);
342  spacing[1] = 2 * M_PI / (sz[1] - 1);
343  Ghost<2, long int> ghost(1);
344 
345  //std::cout << "Spacing: " << spacing[0] << " " << spacing[1] << std::endl;
346 
347  staggered_grid_dist<2, double, aggregate<double, double>> domain(sz, box,ghost,bc);
348 
349  comb<2> c({(char)-1,(char)-1});
350  openfpm::vector<comb<2>> lc;
351  lc.add(c);
352 
353  domain.setStagPosition<phi_>(lc);
354 
355  comb<2> c2({(char)0,(char)-1});
356  lc.clear();
357  lc.add(c2);
358 
359  domain.setStagPosition<ux_>(lc);
360 
361  BOOST_TEST_MESSAGE("Init domain...");
362  auto it = domain.getDomainIterator();
363  while (it.isNext())
364  {
365  auto key_l = it.get();
366  auto key = it.getGKey(key_l);
367  mem_id i = key.get(0);
368  double x = i * spacing[0];
369  mem_id j = key.get(1);
370  double y = j * spacing[1];
371  // Here fill the function value P
372  domain.template getProp<phi_>(key_l) = x*x + 3.0*y*y;
373 
374  ++it;
375  }
376 
377  domain.ghost_get<0>();
378 
379  FD::Derivative_x Dx;
380  FD::Derivative_y Dy;
381 
382  auto ux = FD::getV_stag<ux_>(domain);
383  auto phi = FD::getV_stag<phi_>(domain);
384 
385  ux = Dx(phi) + Dy(phi) + phi + 5;
386 
387  grid_key_dx<2> start({1ul,1u});
388  grid_key_dx<2> stop({(long int)(sz[0]-2),(long int)(sz[1]-2)});
389 
390  double worst_error = 0.0;
391 
392  auto it2 = domain.getSubDomainIterator(start,stop);
393  while (it2.isNext())
394  {
395  auto key = it2.get();
396 
397  double test = domain.template getProp<ux_>(key);
398 
399  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))) );
400  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))) );
401 
402  dx *= 0.5/spacing[0];
403  dy *= 0.5/spacing[1];
404 
405  double phi_inte = 0.5*(domain.template getProp<phi_>(key.move(y,1)) + domain.template getProp<phi_>(key));
406 
407  if (fabs(test - (dx + dy + phi_inte + 5.0)) > worst_error )
408  {worst_error = fabs(test - (dx + dy + phi_inte + 5.0));}
409 
410  ++it2;
411  }
412 
413  BOOST_REQUIRE(worst_error < 5e-13);
414  }
415 
416 
417 BOOST_AUTO_TEST_SUITE_END()
no_equation
Specify the general characteristic of system to solve.
Definition: FD_op_Tests.cpp:46
grid_key_dx
grid_key_dx is the key to access any element in the grid
Definition: grid_key.hpp:18
comb
Position of the element of dimension d in the hyper-cube of dimension dim.
Definition: comb.hpp:34
grid_key_dx::get
__device__ __host__ index_type get(index_type i) const
Get the i index.
Definition: grid_key.hpp:503
Ghost
Definition: Ghost.hpp:39
grid_dist_id
This is a distributed grid.
Definition: grid_dist_id.hpp:96
staggered_grid_dist
Implementation of the staggered grid.
Definition: staggered_dist_grid.hpp:63
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
Point_test::get
auto get() -> decltype(boost::fusion::at_c< i >(data))
getter method for a general property i
Definition: Point_test.hpp:219
Point_test
Test structure used for several test.
Definition: Point_test.hpp:105
openfpm::vector
Implementation of 1-D std::vector like structure.
Definition: map_vector.hpp:202
periodicity
Boundary conditions.
Definition: common.hpp:21