OpenFPM  5.2.0
Project that contain the implementation of distributed structures
DCPSE_op_test_base_tests.cu
1 /*
2  * DCPSE_op_test_base_tests.cpp
3  *
4  * Created on: April 9, 2020
5  * Author: Abhinav Singh
6  *
7  */
8 #include "config.h"
9 #ifdef HAVE_EIGEN
10 #ifdef HAVE_PETSC
11 #define BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS
12 #define BOOST_MPL_LIMIT_VECTOR_SIZE 30
13 
14 
15 
16 #define BOOST_TEST_DYN_LINK
17 
18 #include "util/util_debug.hpp"
19 #include <boost/test/unit_test.hpp>
20 #include <iostream>
21 #include "DCPSE/DCPSE_op/DCPSE_op.hpp"
22 #include "DCPSE/DCPSE_op/DCPSE_Solver.hpp"
23 #include "Operators/Vector/vector_dist_operators.hpp"
24 #include "Vector/vector_dist_subset.hpp"
25 #include "DCPSE/DCPSE_op/EqnsStruct.hpp"
26 #include "DCPSE/DcpseInterpolation.hpp"
27 
28 BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests_cu)
29  BOOST_AUTO_TEST_CASE(dcpse_op_tests) {
30  size_t edgeSemiSize = 40;
31  const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
32  Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
33  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
34  double spacing[2];
35  spacing[0] = 2 * M_PI / (sz[0] - 1);
36  spacing[1] = 2 * M_PI / (sz[1] - 1);
37  Ghost<2, double> ghost(spacing[0] * 3.9);
38  double rCut = 3.9 * spacing[0];
39  BOOST_TEST_MESSAGE("Init vector_dist...");
40  double sigma2 = spacing[0] * spacing[1] / (2 * 4);
41 
42  vector_dist_gpu<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>>> domain(0, box,
43  bc,
44  ghost);
45 
46  //Init_DCPSE(domain)
47  BOOST_TEST_MESSAGE("Init domain...");
48 
49  auto it = domain.getGridIterator(sz);
50  size_t pointId = 0;
51  size_t counter = 0;
52  double minNormOne = 999;
53  while (it.isNext()) {
54  domain.add();
55  auto key = it.get();
56  mem_id k0 = key.get(0);
57  double x = k0 * spacing[0];
58  domain.getLastPos()[0] = x;//+ gaussian(rng);
59  mem_id k1 = key.get(1);
60  double y = k1 * spacing[1];
61  domain.getLastPos()[1] = y;//+gaussian(rng);
62  // Here fill the function value
63  domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
64  domain.template getLastProp<2>() = 2*cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]);
65  ++counter;
66  ++it;
67  }
68  BOOST_TEST_MESSAGE("Sync domain across processors...");
69 
70  domain.map();
71  domain.ghost_get<0>();
72 
73  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
74 
75  Derivative_x_gpu<decltype(verletList)> Dx(domain, verletList, 2, rCut);
76  Derivative_y_gpu<decltype(verletList)> Dy(domain, verletList, 2, rCut);
77  Gradient_gpu<decltype(verletList)> Grad(domain, verletList, 2, rCut);
78  Laplacian_gpu<decltype(verletList)> Lap(domain, verletList, 2, rCut);
79  auto v = getV<1>(domain);
80  auto P = getV<0>(domain);
81 
82  v = 2*Dx(P) + Dy(P);
83  auto it2 = domain.getDomainIterator();
84 
85  double worst = 0.0;
86 
87  while (it2.isNext()) {
88  auto p = it2.get();
89 
90  if (fabs(domain.getProp<1>(p) - domain.getProp<2>(p)) > worst) {
91  worst = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
92  }
93 
94  ++it2;
95  }
96 
97  domain.deleteGhost();
98  BOOST_REQUIRE(worst < 0.03);
99  }
100 
101  BOOST_AUTO_TEST_CASE(dcpse_op_test_lap) {
102  size_t edgeSemiSize = 81;
103  const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
104  Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
105  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
106  double spacing[2];
107  spacing[0] = 2 * M_PI / (sz[0] - 1);
108  spacing[1] = 2 * M_PI / (sz[1] - 1);
109  Ghost<2, double> ghost(spacing[0] * 3.9);
110  double rCut = 3.9 * spacing[0];
111  BOOST_TEST_MESSAGE("Init vector_dist...");
112  double sigma2 = spacing[0] * spacing[1] / (2 * 4);
113 
114  vector_dist_gpu<2, double, aggregate<double, double, double, double, VectorS<2, double>>> domain(0, box,
115  bc,
116  ghost);
117 
118  //Init_DCPSE(domain)
119  BOOST_TEST_MESSAGE("Init domain...");
120  std::normal_distribution<> gaussian{0, sigma2};
121 
122  auto it = domain.getGridIterator(sz);
123  size_t pointId = 0;
124  size_t counter = 0;
125  double minNormOne = 999;
126  while (it.isNext()) {
127  domain.add();
128  auto key = it.get();
129  mem_id k0 = key.get(0);
130  double x = k0 * spacing[0];
131  domain.getLastPos()[0] = x;//+ gaussian(rng);
132  mem_id k1 = key.get(1);
133  double y = k1 * spacing[1];
134  domain.getLastPos()[1] = y;//+gaussian(rng);
135  // Here fill the function value
136  domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
137  // Here fill the validation value for Lap
138  domain.template getLastProp<1>() = - sin(domain.getLastPos()[0]) - sin(domain.getLastPos()[1]);
139  domain.template getLastProp<4>()[0] = -sin(domain.getLastPos()[0]);
140  domain.template getLastProp<4>()[1] = -sin(domain.getLastPos()[1]);
141 
142  ++counter;
143  ++it;
144  }
145  BOOST_TEST_MESSAGE("Sync domain across processors...");
146 
147  domain.map();
148  domain.ghost_get<0>();
149 
150  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
151 
152  Laplacian_gpu<decltype(verletList)> Lap(domain, verletList, 2, rCut);
153  auto v = getV<1>(domain);
154  auto P = getV<0>(domain);
155  auto vv = getV<2>(domain);
156  auto errv = getV<3>(domain);
157 
158  vv = Lap(P);
159  auto it2 = domain.getDomainIterator();
160 
161  double worst = 0.0;
162 
163  while (it2.isNext()) {
164  auto p = it2.get();
165 
166  if (fabs(domain.getProp<1>(p) - domain.getProp<2>(p)) > worst) {
167  worst = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
168  }
169 
170  ++it2;
171  }
172  errv=v-vv;
173  domain.deleteGhost();
174  BOOST_REQUIRE(worst < 0.3);
175  }
176 
177  BOOST_AUTO_TEST_CASE(dcpse_op_div) {
178 // int rank;
179 // MPI_Comm_rank(MPI_COMM_WORLD, &rank);
180  size_t edgeSemiSize = 31;
181  const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
182  Box<2, double> box({0, 0}, {10.0, 10.0});
183  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
184  double spacing[2];
185  spacing[0] = box.getHigh(0) / (sz[0] - 1);
186  spacing[1] = box.getHigh(1) / (sz[1] - 1);
187  double rCut = 3.1* spacing[0];
188  Ghost<2, double> ghost(rCut);
189  BOOST_TEST_MESSAGE("Init vector_dist...");
190  double sigma2 = spacing[0] * spacing[1] / (2 * 4);
191 
192  vector_dist_gpu<2, double, aggregate<double, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>,double,double>> domain(
193  0, box, bc, ghost);
194 
195  //Init_DCPSE(domain)
196  BOOST_TEST_MESSAGE("Init domain...");
197 // std::random_device rd{};
198 // std::mt19937 rng{rd()};
199  std::mt19937 rng{6666666};
200 
201  std::normal_distribution<> gaussian{0, sigma2};
202 
203  auto it = domain.getGridIterator(sz);
204  size_t pointId = 0;
205  size_t counter = 0;
206  double minNormOne = 999;
207  while (it.isNext()) {
208  domain.add();
209  auto key = it.get();
210  mem_id k0 = key.get(0);
211  double x = k0 * spacing[0];
212  domain.getLastPos()[0] = x;//+ gaussian(rng);
213  mem_id k1 = key.get(1);
214  double y = k1 * spacing[1];
215  domain.getLastPos()[1] = y;//+gaussian(rng);
216  // Here fill the function value
217  domain.template getLastProp<1>()[0] = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
218  domain.template getLastProp<1>()[1] = cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]);
219 
220 
221  // Here fill the validation value for Divergence
222  domain.template getLastProp<0>()= cos(domain.getLastPos()[0]) - sin(domain.getLastPos()[1]);
223  /* domain.template getLastProp<4>()[0] =
224  cos(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) +
225  cos(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
226  domain.template getLastProp<4>()[1] =
227  -sin(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) -
228  sin(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));*/
229 
230 
231  ++counter;
232  ++it;
233  }
234  BOOST_TEST_MESSAGE("Sync domain across processors...");
235 
236  domain.map();
237  domain.ghost_get<0>();
238 
239  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
240 
241  Divergence_gpu<decltype(verletList)> Div(domain, verletList, 2, rCut);
242  Derivative_x_gpu<decltype(verletList)> Dx(domain, verletList, 2, rCut);
243  Derivative_y_gpu<decltype(verletList)> Dy(domain, verletList, 2, rCut);
244 
245  auto v = getV<1>(domain);
246  auto anasol = getV<0>(domain);
247  auto div = getV<5>(domain);
248  auto div2 = getV<6>(domain);
249 
250  domain.ghost_get<1>();
251  div = Div(v);
252  // div2=Dx(v[0])+Dy(v[1]);
253  // auto it2 = domain.getDomainIterator();
254 
255  // double worst1 = 0.0;
256  // double worst2 = 0.0;
257 
258  // while (it2.isNext()) {
259  // auto p = it2.get();
260  // if (fabs(domain.getProp<0>(p) - domain.getProp<5>(p)) > worst1) {
261  // worst1 = fabs(domain.getProp<0>(p) - domain.getProp<5>(p));
262  // }
263  // if (fabs(domain.getProp<0>(p) - domain.getProp<6>(p)) > worst2) {
264  // worst2 = fabs(domain.getProp<0>(p) - domain.getProp<6>(p));
265  // }
266  // ++it2;
267  // }
268 
269  domain.deleteGhost();
270 /*
271  domain.write("DIV");
272 
273  std::cout<<worst1<<":"<<worst2;
274 */
275 
276  // BOOST_REQUIRE(worst1 < 0.05);
277  // BOOST_REQUIRE(worst2 < 0.05);
278 
279  }
280 
281  BOOST_AUTO_TEST_CASE(dcpse_op_vec) {
282  // int rank;
283  // MPI_Comm_rank(MPI_COMM_WORLD, &rank);
284  size_t edgeSemiSize = 160;
285  const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
286  Box<2, double> box({0, 0}, {10,10});
287  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
288  double spacing[2];
289  spacing[0] = box.getHigh(0) / (sz[0] - 1);
290  spacing[1] = box.getHigh(1) / (sz[1] - 1);
291  Ghost<2, double> ghost(spacing[0] * 3.1);
292  double rCut = 3.1 * spacing[0];
293  BOOST_TEST_MESSAGE("Init vector_dist...");
294  double sigma2 = spacing[0] * spacing[1] / (2 * 4);
295 
296  vector_dist_gpu<2, double, aggregate<double, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>,double,double>> domain(
297  0, box, bc, ghost);
298 
299  //Init_DCPSE(domain)
300  BOOST_TEST_MESSAGE("Init domain...");
301  std::random_device rd{};
302  // std::mt19937 rng{rd()};
303  std::mt19937 rng{6666666};
304 
305  std::normal_distribution<> gaussian{0, sigma2};
306 
307  auto it = domain.getGridIterator(sz);
308  size_t pointId = 0;
309  size_t counter = 0;
310  double minNormOne = 999;
311  while (it.isNext()) {
312  domain.add();
313  auto key = it.get();
314  mem_id k0 = key.get(0);
315  double x = k0 * spacing[0];
316  domain.getLastPos()[0] = x;//+ gaussian(rng);
317  mem_id k1 = key.get(1);
318  double y = k1 * spacing[1];
319  domain.getLastPos()[1] = y;//+gaussian(rng);
320  // Here fill the function value
321  domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
322  domain.template getLastProp<1>()[0] = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
323  domain.template getLastProp<1>()[1] = cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]);
324 
325  // Here fill the validation value for Df/Dx
326  domain.template getLastProp<2>()[0] = 0;
327  domain.template getLastProp<2>()[1] = 0;
328  domain.template getLastProp<4>()[0] =
329  cos(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) +
330  cos(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
331  domain.template getLastProp<4>()[1] =
332  -sin(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) -
333  sin(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
334 
335  domain.template getLastProp<5>() = cos(domain.getLastPos()[0])*(sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]))+cos(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
336 
337 
338  ++counter;
339  ++it;
340  }
341  BOOST_TEST_MESSAGE("Sync domain across processors...");
342 
343  domain.map();
344  domain.ghost_get<0>();
345 
346  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
347 
348  Advection_gpu<decltype(verletList)> Adv(domain, verletList, 2, rCut);
349  auto v = getV<1>(domain);
350  auto P = getV<0>(domain);
351  auto dv = getV<3>(domain);
352  auto dP = getV<6>(domain);
353 
354 
355  domain.ghost_get<1>();
356  dv = Adv(v, v);
357  auto it2 = domain.getDomainIterator();
358 
359  double worst1 = 0.0;
360 
361  while (it2.isNext()) {
362  auto p = it2.get();
363 
364  if (fabs(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]) > worst1) {
365  worst1 = fabs(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]);
366 
367  }
368 
369  ++it2;
370  }
371 
372  //std::cout << "Maximum Error in component 2: " << worst1 << std::endl;
373  //domain.write("v1");
374  BOOST_REQUIRE(worst1 < 0.03);
375 
376 
377  dP = Adv(v, P);
378  auto it3 = domain.getDomainIterator();
379 
380  double worst2 = 0.0;
381 
382  while (it3.isNext()) {
383  auto p = it3.get();
384  if (fabs(domain.getProp<6>(p) - domain.getProp<5>(p)) > worst2) {
385  worst2 = fabs(domain.getProp<6>(p) - domain.getProp<5>(p));
386 
387  }
388 
389  ++it3;
390  }
391 
392 
393  domain.deleteGhost();
394  //domain.write("v2");
395  BOOST_REQUIRE(worst2 < 0.03);
396 
397  }
398 
399 
400  BOOST_AUTO_TEST_CASE(dcpse_slice) {
401  const size_t sz[2] = {31,31};
402  Box<2, double> box({0, 0}, {1,1});
403  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
404  double spacing = box.getHigh(0) / (sz[0] - 1);
405  Ghost<2, double> ghost(spacing * 3.9);
406  double rCut = 3.9 * spacing;
407 
408  vector_dist_gpu<2, double, aggregate<double,VectorS<2, double>,double,double[3][3]>> Particles(0, box, bc, ghost);
409 
410  auto it = Particles.getGridIterator(sz);
411  while (it.isNext()) {
412  Particles.add();
413  auto key = it.get();
414  double x = key.get(0) * it.getSpacing(0);
415  Particles.getLastPos()[0] = x;
416  double y = key.get(1) * it.getSpacing(1);
417  Particles.getLastPos()[1] = y;
418 
419  Particles.getLastProp<1>()[0] = sin(x+y);
420  Particles.getLastProp<1>()[1] = cos(x+y);
421 
422  ++it;
423  }
424 
425  Particles.map();
426  Particles.ghost_get<0,1>();
427 
428  auto P = getV<0>(Particles);
429  auto V = getV<1>(Particles);
430  auto S = getV<2>(Particles);
431  auto Sig = getV<3>(Particles);
432 
433  auto verletList = Particles.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
434 
435  Derivative_x_gpu<decltype(verletList)> Dx(Particles, verletList, 2, rCut);
436 
437  P = Dx(V[0]);
438  S = V[0]*V[0] + V[1]*V[1];
439 
440  Sig[0][1] = V[0]*V[0] + V[1]*V[1];
441  Sig[1][0] = P;
442  Sig[2][2] = 5.0;
443 
444  auto it2 = Particles.getDomainIterator();
445 
446  double err = 0.0;
447 
448  while (it2.isNext())
449  {
450  auto p = it2.get();
451 
452  if (fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]) >= err )
453  {
454  err = fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]);
455  }
456 
457  if (fabs(Particles.getProp<2>(p) - 1.0) >= err )
458  {
459  err = fabs(Particles.getProp<2>(p) - 1.0);
460  }
461 
462  if (fabs(Particles.getProp<3>(p)[0][1] - 1.0) >= err )
463  {
464  err = fabs(Particles.getProp<3>(p)[0][1] - 1.0);
465  }
466 
467  if (fabs(Particles.getProp<3>(p)[1][0] - Particles.getProp<1>(p)[1]) >= err )
468  {
469  err = fabs(Particles.getProp<3>(p)[1][0] - Particles.getProp<1>(p)[1]);
470  }
471 
472  if (fabs(Particles.getProp<3>(p)[2][2] - 5.0) >= err )
473  {
474  err = fabs(Particles.getProp<3>(p)[2][2] - 5.0);
475  }
476 
477  ++it2;
478  }
479 
480  //Particles.write("test_out");
481  //std::cout << "Error: " << err << " " << create_vcluster().rank() << std::endl;
482  BOOST_REQUIRE(err < 0.03);
483 }
484 
485  BOOST_AUTO_TEST_CASE(dcpse_slice_3d) {
486  const size_t sz[3] = {17,17,17};
487  Box<3, double> box({0, 0,0}, {1,1,1});
488  size_t bc[3] = {NON_PERIODIC, NON_PERIODIC,NON_PERIODIC};
489  double spacing = box.getHigh(0) / (sz[0] - 1);
490  Ghost<3, double> ghost(spacing * 3.9);
491  double rCut = 3.9 * spacing;
492 
493  vector_dist_gpu<3, double, aggregate<double,VectorS<3, double>,double,double[3][3]>> Particles(0, box, bc, ghost);
494 
495  auto it = Particles.getGridIterator(sz);
496  while (it.isNext()) {
497  Particles.add();
498  auto key = it.get();
499  double x = key.get(0) * it.getSpacing(0);
500  Particles.getLastPos()[0] = x;
501  double y = key.get(1) * it.getSpacing(1);
502  Particles.getLastPos()[1] = y;
503  double z = key.get(2) * it.getSpacing(2);
504  Particles.getLastPos()[2] = z;
505 
506  Particles.getLastProp<1>()[0] = sin(x+y);
507  Particles.getLastProp<1>()[1] = cos(x+y);
508  Particles.getLastProp<1>()[2] = 1.0;
509 
510  ++it;
511  }
512 
513  Particles.map();
514  Particles.ghost_get<0,1>();
515 
516 
517  auto P = getV<0>(Particles);
518  auto V = getV<1>(Particles);
519  auto S = getV<2>(Particles);
520  auto Sig = getV<3>(Particles);
521 
522  auto verletList = Particles.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
523 
524  Derivative_x_gpu<decltype(verletList)> Dx(Particles, verletList, 2, rCut);
525 
526  P = Dx(V[0]);
527  S = V[0]*V[0] + V[1]*V[1]+V[2]*V[2];
528 
529  Sig[0][1] = V[0]*V[0] + V[1]*V[1]+V[2]*V[2];
530  Sig[1][0] = P;
531  Sig[2][2] = 5.0;
532 
533  auto it2 = Particles.getDomainIterator();
534 
535  double err1 = 0.0;
536  double err2 = 0.0;
537  double err3 = 0.0;
538  double err4 = 0.0;
539  double err5 = 0.0;
540 
541  while (it2.isNext())
542  {
543  auto p = it2.get();
544 
545  // Here we check that P = Dx(V[0]) works P is Prop 0 = sin(x+y) V[0] = sin(x+y) and V[1] is cos(x+y)
546  if (fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]) >= err1 )
547  {
548  err1 = fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]);
549  }
550 
551  if (fabs(Particles.getProp<2>(p) - 2.0) >= err2 )
552  {
553  err2 = fabs(Particles.getProp<2>(p) - 2.0);
554  }
555 
556  if (fabs(Particles.getProp<3>(p)[0][1] - 2.0) >= err3 )
557  {
558  err3 = fabs(Particles.getProp<3>(p)[0][1] - 2.0);
559  }
560 
561  // V[0]*V[0] + V[1]*V[1]+V[2]*V[2] = 2 and
562  if (fabs(Particles.getProp<3>(p)[0][1] - Particles.getProp<2>(p)) >= err4 )
563  {
564  err4 = fabs(Particles.getProp<3>(p)[0][1] - Particles.getProp<2>(p));
565  }
566 
567  if (fabs(Particles.getProp<3>(p)[2][2] - 5.0) >= err5 )
568  {
569  err5 = fabs(Particles.getProp<3>(p)[2][2] - 5.0);
570  }
571 
572  ++it2;
573  }
574 
575  //std::cout << err1 << " " << err2 << " " << err3 << " " << err4 << " " << err5 << std::endl;
576  BOOST_REQUIRE(err1 < 0.08);
577  BOOST_REQUIRE(err2 < 0.03);
578  BOOST_REQUIRE(err3 < 0.03);
579  BOOST_REQUIRE(err4 < 0.03);
580  BOOST_REQUIRE(err5 < 0.03);
581  }
582 
583  BOOST_AUTO_TEST_CASE(dcpse_op_convection) {
584  size_t edgeSemiSize = 20;
585  const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
586  Box<2, double> box({-1, -1}, {1.0, 1.0});
587  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
588  double spacing[2];
589  spacing[0] = 2.0/ (sz[0] - 1);
590  spacing[1] = 2.0 / (sz[1] - 1);
591  double rCut = 3.1 * spacing[0];
592  Ghost<2, double> ghost(rCut);
593 
594  BOOST_TEST_MESSAGE("Init vector_dist...");
595 
596  vector_dist_gpu<2, double, aggregate<double, double, double, VectorS<2, double>>> domain(0, box,bc,
597  ghost);
598  domain.setPropNames({"Concentration","Concentration_temp","Temp","Velocity"});
599 
600  //Init_DCPSE(domain)
601  BOOST_TEST_MESSAGE("Init domain...");
602 
603  auto it = domain.getGridIterator(sz);
604  while (it.isNext()) {
605  domain.add();
606  auto key = it.get();
607  mem_id k0 = key.get(0);
608  double x = -1.0+k0 * spacing[0];
609  domain.getLastPos()[0] = x;//+ gaussian(rng);
610  mem_id k1 = key.get(1);
611  double y = -1.0+k1 * spacing[1];
612  domain.getLastPos()[1] = y;//+gaussian(rng);
613  // Here fill the function value
614  if (x>-1 && y>-1 && x<1 && y<1)
615  {
616  domain.template getLastProp<3>()[0] = (-y)*exp(-10*((x)*(x)+(y)*(y)));;
617  domain.template getLastProp<3>()[1] = (x)*exp(-10*((x)*(x)+(y)*(y)));;
618  }
619  else{
620  domain.template getLastProp<3>()[0] = 0.0;
621  domain.template getLastProp<3>()[1] = 0.0;
622  }
623  if (x==0.0 && y>-0.5 && y<0.5)
624  {
625  domain.template getLastProp<0>() = 1.0;
626  }
627  else
628  {
629  domain.template getLastProp<0>() = 0.0;
630  }
631  ++it;
632  }
633  BOOST_TEST_MESSAGE("Sync domain across processors...");
634 
635  domain.map();
636  domain.ghost_get<0>();
637 
638  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
639 
640  //Derivative_x_gpu Dx(domain, 2, rCut);
641  Derivative_xx_gpu<decltype(verletList)> Dxx(domain, verletList, 2, rCut);
642  //Derivative_y_gpu Dy(domain, 2, rCut);
643  Derivative_yy_gpu<decltype(verletList)> Dyy(domain, verletList, 2, rCut);
644  auto C = getV<0>(domain);
645  auto V = getV<3>(domain);
646  auto Cnew = getV<1>(domain);
647  auto Pos = getV<POS_PROP>(domain);
648  timer tt;
649  //domain.write_frame("Convection_init",0);
650  int ctr=0;
651  double t=0,tf=1,dt=1e-2;
652  while(t<tf)
653  {
654  domain.write_frame("Convection",ctr);
655  domain.ghost_get<0>();
656  Cnew=C+dt*0.01*(Dxx(C)+Dyy(C));
657  C=Cnew;
658  Pos=Pos+dt*V;
659  domain.map();
660  domain.ghost_get<0>();
661  auto it2 = domain.getDomainIterator();
662  while (it2.isNext()) {
663  auto p = it2.get();
664  Point<2, double> xp = domain.getPos(p);
665  double x=xp[0],y=xp[1];
666  if (x>-1 && y>-1 && x<1 && y<1)
667  {
668  domain.template getProp<3>(p)[0] = (-y)*exp(-10*((x)*(x)+(y)*(y)));
669  domain.template getProp<3>(p)[1] = (x)*exp(-10*((x)*(x)+(y)*(y)));;
670  }
671  else{
672  domain.template getProp<3>(p)[0] = 0.0;
673  domain.template getProp<3>(p)[1] = 0.0;
674  }
675  ++it2;
676  }
677  tt.start();
678  Dxx.update(domain);
679  Dyy.update(domain);
680  tt.stop();
681  //std::cout<tt.getwct()<<""
682  ctr++;
683  t+=dt;
684  }
685 
686  Dxx.deallocate(domain);
687  Dyy.deallocate(domain);
688 
689  /* std::cout<<"Dx"<<std::endl;
690  Dx.checkMomenta(domain);
691  std::cout<<"Dy"<<std::endl;
692  Dy.checkMomenta(domain);
693  std::cout<<"Dxx"<<std::endl;
694  Dxx.checkMomenta(domain);
695  int ctr=0;
696  P=0;
697  v=0;
698  K1=0;
699  auto its2 = domain.getDomainIterator();
700  while (its2.isNext()) {
701  auto p = its2.get();
702  Dx.DrawKernel<0>(domain, p.getKey());
703  Dy.DrawKernel<1>(domain, p.getKey());
704  Dxx.DrawKernel<2>(domain, p.getKey());
705  domain.write_frame("Kernels",ctr);
706  P=0;
707  v=0;
708  K1=0;
709  ++its2;
710  ctr++;
711  }*/
712 
713  }
714 
715 
716 
717 BOOST_AUTO_TEST_SUITE_END()
718 
719 
720 #endif
721 #endif
722 
723 
724 
Box
This class represent an N-dimensional box.
Definition: Box.hpp:60
Ghost
Definition: Ghost.hpp:40
Lap
Laplacian second order on h (spacing)
Definition: Laplacian.hpp:23
Point_test
Test structure used for several test.
Definition: Point_test.hpp:106
Point
This class implement the point shape in an N-dimensional space.
Definition: Point.hpp:28
timer
Class for cpu time benchmarking.
Definition: timer.hpp:28
timer::stop
void stop()
Stop the timer.
Definition: timer.hpp:119
timer::start
void start()
Start the timer.
Definition: timer.hpp:90
vector_dist
Distributed vector.
Definition: vector_dist.hpp:176