OpenFPM  5.2.0
Project that contain the implementation of distributed structures
DCPSE_op_test_base_tests.cpp
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 "Vector/vector_dist_multiphase_functions.hpp"
26 #include "DCPSE/DCPSE_op/EqnsStruct.hpp"
27 #include "DCPSE/DcpseInterpolation.hpp"
28 
29 BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
30 BOOST_AUTO_TEST_CASE(dcpse_op_tests) {
31  size_t edgeSemiSize = 40;
32  const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
33  Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
34  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
35  double spacing[2];
36  spacing[0] = 2 * M_PI / (sz[0] - 1);
37  spacing[1] = 2 * M_PI / (sz[1] - 1);
38  Ghost<2, double> ghost(spacing[0] * 3.9);
39  double rCut = 3.9 * spacing[0];
40  BOOST_TEST_MESSAGE("Init vector_dist...");
41  double sigma2 = spacing[0] * spacing[1] / (2 * 4);
42 
43  vector_dist<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>>> domain(0, box,
44  bc,
45  ghost);
46 
47  //Init_DCPSE(domain)
48  BOOST_TEST_MESSAGE("Init domain...");
49 
50  auto it = domain.getGridIterator(sz);
51  size_t pointId = 0;
52  size_t counter = 0;
53  double minNormOne = 999;
54  while (it.isNext()) {
55  domain.add();
56  auto key = it.get();
57  mem_id k0 = key.get(0);
58  double x = k0 * spacing[0];
59  domain.getLastPos()[0] = x;//+ gaussian(rng);
60  mem_id k1 = key.get(1);
61  double y = k1 * spacing[1];
62  domain.getLastPos()[1] = y;//+gaussian(rng);
63  // Here fill the function value
64  domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
65  domain.template getLastProp<2>() = 2*cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]);
66  ++counter;
67  ++it;
68  }
69  BOOST_TEST_MESSAGE("Sync domain across processors...");
70 
71  domain.map();
72  domain.ghost_get<0>();
73 
74  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
75 
76  Derivative_x<decltype(verletList)> Dx(domain, verletList, 2, rCut);
77  Derivative_y<decltype(verletList)> Dy(domain, verletList, 2, rCut);
78  Gradient<decltype(verletList)> Grad(domain, verletList, 2, rCut);
79  Laplacian<decltype(verletList)> Lap(domain, verletList, 2, rCut);
80  auto v = getV<1>(domain);
81  auto P = getV<0>(domain);
82 
83  v = 2*Dx(P) + Dy(P);
84  auto it2 = domain.getDomainIterator();
85 
86  double worst = 0.0;
87 
88  while (it2.isNext()) {
89  auto p = it2.get();
90 
91  if (fabs(domain.getProp<1>(p) - domain.getProp<2>(p)) > worst) {
92  worst = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
93  }
94 
95  ++it2;
96  }
97 
98  domain.deleteGhost();
99  BOOST_REQUIRE(worst < 0.03);
100 
101  }
102 
103  BOOST_AUTO_TEST_CASE(dcpse_op_save_load) {
104  size_t edgeSemiSize = 40;
105  const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
106  Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
107  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
108  double spacing[2];
109  spacing[0] = 2 * M_PI / (sz[0] - 1);
110  spacing[1] = 2 * M_PI / (sz[1] - 1);
111  Ghost<2, double> ghost(spacing[0] * 3.9);
112  double rCut = 3.9 * spacing[0];
113  BOOST_TEST_MESSAGE("Init vector_dist...");
114  double sigma2 = spacing[0] * spacing[1] / (2 * 4);
115 
116  vector_dist<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>>> domain(0, box,
117  bc,
118  ghost);
119 
120  //Init_DCPSE(domain)
121  BOOST_TEST_MESSAGE("Init domain...");
122 
123  auto it = domain.getGridIterator(sz);
124  size_t pointId = 0;
125  size_t counter = 0;
126  double minNormOne = 999;
127  while (it.isNext()) {
128  domain.add();
129  auto key = it.get();
130  mem_id k0 = key.get(0);
131  double x = k0 * spacing[0];
132  domain.getLastPos()[0] = x;//+ gaussian(rng);
133  mem_id k1 = key.get(1);
134  double y = k1 * spacing[1];
135  domain.getLastPos()[1] = y;//+gaussian(rng);
136  // Here fill the function value
137  domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
138  domain.template getLastProp<2>() = 2*cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]);
139  ++counter;
140  ++it;
141  }
142  BOOST_TEST_MESSAGE("Sync domain across processors...");
143 
144  domain.map();
145  domain.ghost_get<0>();
146 
147  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
148  Derivative_x<decltype(verletList)> Dx(domain, verletList, 2, rCut);
149  Derivative_y<decltype(verletList)> Dy(domain, verletList, 2, rCut);
150  auto v = getV<1>(domain);
151  auto v2 = getV<3>(domain);
152  auto P = getV<0>(domain);
153  v2 = 2*Dx(P) + Dy(P);
154  Dx.save(domain,"DX_test");
155  Dy.save(domain,"DY_test");
156  Derivative_x<decltype(verletList)> DxLoaded(domain, verletList, 2, rCut, support_options::LOAD);
157  Derivative_y<decltype(verletList)> DyLoaded(domain, verletList, 2, rCut, support_options::LOAD);
158  DxLoaded.load(domain,"DX_test");
159  DyLoaded.load(domain,"DY_test");
160  v= 2*DxLoaded(P)+DyLoaded(P);
161  auto it2 = domain.getDomainIterator();
162  double worst = 0.0;
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  domain.deleteGhost();
173  //std::cout<<worst;
174  BOOST_REQUIRE(worst < 0.03);
175  }
176 
177  BOOST_AUTO_TEST_CASE(dcpse_op_save_load2) {
178  size_t edgeSemiSize = 40;
179  const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
180  Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
181  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
182  double spacing[2];
183  spacing[0] = 2 * M_PI / (sz[0] - 1);
184  spacing[1] = 2 * M_PI / (sz[1] - 1);
185  Ghost<2, double> ghost(spacing[0] * 3.9);
186  double rCut = 3.9 * spacing[0];
187  BOOST_TEST_MESSAGE("Init vector_dist...");
188  double sigma2 = spacing[0] * spacing[1] / (2 * 4);
189 
190  vector_dist<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>>> domain(0, box,
191  bc,
192  ghost);
193 
194  //Init_DCPSE(domain)
195  BOOST_TEST_MESSAGE("Init domain...");
196 
197  auto it = domain.getGridIterator(sz);
198  size_t pointId = 0;
199  size_t counter = 0;
200  double minNormOne = 999;
201  while (it.isNext()) {
202  domain.add();
203  auto key = it.get();
204  mem_id k0 = key.get(0);
205  double x = k0 * spacing[0];
206  domain.getLastPos()[0] = x;//+ gaussian(rng);
207  mem_id k1 = key.get(1);
208  double y = k1 * spacing[1];
209  domain.getLastPos()[1] = y;//+gaussian(rng);
210  // Here fill the function value
211  domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
212  domain.template getLastProp<2>() = 2*cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]);
213  ++counter;
214  ++it;
215  }
216  BOOST_TEST_MESSAGE("Sync domain across processors...");
217 
218  domain.map();
219  domain.ghost_get<0>();
220  auto v = getV<1>(domain);
221  auto v2 = getV<3>(domain);
222  auto P = getV<0>(domain);
223  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
224  Derivative_x<decltype(verletList)> DxLoaded(domain, verletList, 2, rCut, support_options::LOAD);
225  Derivative_y<decltype(verletList)> DyLoaded(domain, verletList, 2, rCut, support_options::LOAD);
226  DxLoaded.load(domain,"DX_test");
227  DyLoaded.load(domain,"DY_test");
228  v= 2*DxLoaded(P)+DyLoaded(P);
229  auto it2 = domain.getDomainIterator();
230  double worst = 0.0;
231  while (it2.isNext()) {
232  auto p = it2.get();
233 
234  if (fabs(domain.getProp<1>(p) - domain.getProp<2>(p)) > worst) {
235  worst = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
236  }
237 
238  ++it2;
239  }
240  domain.deleteGhost();
241  //std::cout<<worst;
242  BOOST_REQUIRE(worst < 0.03);
243  }
244 
245  BOOST_AUTO_TEST_CASE(dcpse_op_tests_fa) {
246  size_t edgeSemiSize = 40;
247  const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
248  Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
249  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
250  double spacing[2];
251  spacing[0] = 2 * M_PI / (sz[0] - 1);
252  spacing[1] = 2 * M_PI / (sz[1] - 1);
253  Ghost<2, double> ghost(spacing[0] * 3.9);
254  double rCut = 4.1 * spacing[0];
255  BOOST_TEST_MESSAGE("Init vector_dist...");
256  double sigma2 = spacing[0] * spacing[1] / (2 * 4);
257 
258  typedef vector_dist<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>>> vector_type;
259 
260  vector_type domain(0, box,bc,ghost);
261 
262  //Init_DCPSE(domain)
263  BOOST_TEST_MESSAGE("Init domain...");
264 
265  auto it = domain.getGridIterator(sz);
266  size_t pointId = 0;
267  size_t counter = 0;
268  double minNormOne = 999;
269  while (it.isNext()) {
270  domain.add();
271  auto key = it.get();
272  mem_id k0 = key.get(0);
273  double x = k0 * spacing[0];
274  domain.getLastPos()[0] = x;//+ gaussian(rng);
275  mem_id k1 = key.get(1);
276  double y = k1 * spacing[1];
277  domain.getLastPos()[1] = y;//+gaussian(rng);
278  // Here fill the function value
279  domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
280  domain.template getLastProp<2>() = cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]);
281  ++counter;
282  ++it;
283  }
284  BOOST_TEST_MESSAGE("Sync domain across processors...");
285 
286  domain.map();
287  domain.ghost_get<0>();
288 
289  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC>(rCut);
290  PPInterpolation<vector_type,vector_type,decltype(verletList)> Fx(domain, domain, verletList, 2, rCut);
291  auto v = getV<1>(domain);
292  auto P = getV<0>(domain);
293 
294  Fx.p2p<0,1>();
295  auto it2 = domain.getDomainIterator();
296  double worst = 0.0;
297  while (it2.isNext()) {
298  auto p = it2.get();
299  if (fabs(domain.getProp<1>(p) - domain.getProp<0>(p)) > worst) {
300  worst = fabs(domain.getProp<1>(p) - domain.getProp<0>(p));
301  }
302  ++it2;
303  }
304  //std::cout<<"Worst:"<<worst<<std::endl;
305  domain.deleteGhost();
306  //domain.write_frame("test",0,0.024,BINARY);
307  BOOST_REQUIRE(worst < 0.03);
308  }
309 
310  BOOST_AUTO_TEST_CASE(dcpse_op_tests_mfa) {
311  size_t edgeSemiSize = 40;
312  const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
313  Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
314  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
315  double spacing[2];
316  spacing[0] = 2 * M_PI / (sz[0] - 1);
317  spacing[1] = 2 * M_PI / (sz[1] - 1);
318  Ghost<2, double> ghost(spacing[0] * 3.9);
319  double rCut = 3.9 * spacing[0];
320  BOOST_TEST_MESSAGE("Init vector_dist...");
321  double sigma2 = spacing[0] * spacing[1] / ( 4);
322  std::normal_distribution<> gaussian{0, sigma2};
323  std::mt19937 rng{6666666};
324  typedef vector_dist<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>>> vector_dist1;
325  typedef vector_dist<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>>> vector_dist2;
326 
327  vector_dist1 domain(0, box,bc,ghost);
328  vector_dist2 domain2(domain.getDecomposition(),0);
329 
330  //Init_DCPSE(domain)
331  BOOST_TEST_MESSAGE("Init domain...");
332 
333  auto it = domain.getGridIterator(sz);
334  size_t pointId = 0;
335  size_t counter = 0;
336  double minNormOne = 999;
337  while (it.isNext()) {
338  domain.add();
339  domain2.add();
340  auto key = it.get();
341  mem_id k0 = key.get(0);
342  mem_id k1 = key.get(1);
343  double x = k0 * spacing[0];
344  double y = k1 * spacing[1];
345  domain.getLastPos()[0] = x;//+ gaussian(rng);
346  domain.getLastPos()[1] = y;//+gaussian(rng);
347  if(x!=0 && y!=0 && x!=box.getHigh(0) && y!=box.getHigh(1)){
348  domain2.getLastPos()[0] = x+ gaussian(rng);
349  domain2.getLastPos()[1] = y+ gaussian(rng);
350  }
351  else{
352  domain2.getLastPos()[0] = x;
353  domain2.getLastPos()[1] = y;
354  }
355  // Here fill the function value
356  domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
357  domain.template getLastProp<1>() = 0.0;
358  domain2.template getLastProp<0>() = sin(domain2.getLastPos()[0]) + sin(domain2.getLastPos()[1]);
359  ++counter;
360  ++it;
361  }
362  BOOST_TEST_MESSAGE("Sync domain across processors...");
363 
364  domain.map();
365  domain2.map();
366  domain.ghost_get<0>();
367  domain2.ghost_get<0>();
368 
369  auto cellListDomain2 = domain2.getCellList(rCut);
370  auto verletList = createVerletTwoPhase(domain,domain2,cellListDomain2,rCut);
371 
372  PPInterpolation<vector_dist2,vector_dist1,decltype(verletList)> Fx(domain2, domain, verletList, 2, rCut);
373  //auto v = getV<1>(domain);
374  //auto P = getV<0>(domain);
375  Fx.p2p<0,1>();
376  auto it2 = domain.getDomainIterator();
377  double worst = 0.0;
378  while (it2.isNext()) {
379  auto p = it2.get();
380  //domain.template getProp<2>(p) = domain.getProp<1>(p) - domain.getProp<0>(p);
381  if (fabs(domain.getProp<1>(p) - domain.getProp<0>(p)) > worst) {
382  worst = fabs(domain.getProp<1>(p) - domain.getProp<0>(p));
383  }
384  ++it2;
385  }
386  //std::cout<<"Worst:"<<worst<<std::endl;
387  domain.deleteGhost();
388  //domain.write("test1");
389  //domain2.write("test2");
390  BOOST_REQUIRE(worst < 0.03);
391  }
392 
393 
394  BOOST_AUTO_TEST_CASE(dcpse_op_test_lap) {
395  size_t edgeSemiSize = 81;
396  const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
397  Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
398  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
399  double spacing[2];
400  spacing[0] = 2 * M_PI / (sz[0] - 1);
401  spacing[1] = 2 * M_PI / (sz[1] - 1);
402  Ghost<2, double> ghost(spacing[0] * 3.9);
403  double rCut = 3.9 * spacing[0];
404  BOOST_TEST_MESSAGE("Init vector_dist...");
405  double sigma2 = spacing[0] * spacing[1] / (2 * 4);
406 
407  vector_dist<2, double, aggregate<double, double, double, double, VectorS<2, double>>> domain(0, box,
408  bc,
409  ghost);
410 
411  //Init_DCPSE(domain)
412  BOOST_TEST_MESSAGE("Init domain...");
413  std::normal_distribution<> gaussian{0, sigma2};
414 
415  auto it = domain.getGridIterator(sz);
416  size_t pointId = 0;
417  size_t counter = 0;
418  double minNormOne = 999;
419  while (it.isNext()) {
420  domain.add();
421  auto key = it.get();
422  mem_id k0 = key.get(0);
423  double x = k0 * spacing[0];
424  domain.getLastPos()[0] = x;//+ gaussian(rng);
425  mem_id k1 = key.get(1);
426  double y = k1 * spacing[1];
427  domain.getLastPos()[1] = y;//+gaussian(rng);
428  // Here fill the function value
429  domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
430  // Here fill the validation value for Lap
431  domain.template getLastProp<1>() = - sin(domain.getLastPos()[0]) - sin(domain.getLastPos()[1]);
432  domain.template getLastProp<4>()[0] = -sin(domain.getLastPos()[0]);
433  domain.template getLastProp<4>()[1] = -sin(domain.getLastPos()[1]);
434 
435  ++counter;
436  ++it;
437  }
438  BOOST_TEST_MESSAGE("Sync domain across processors...");
439 
440  domain.map();
441  domain.ghost_get<0>();
442 
443  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
444  Laplacian<decltype(verletList)> Lap(domain, verletList, 2, rCut);
445  auto v = getV<1>(domain);
446  auto P = getV<0>(domain);
447  auto vv = getV<2>(domain);
448  auto errv = getV<3>(domain);
449 
450  vv = Lap(P);
451  auto it2 = domain.getDomainIterator();
452 
453  double worst = 0.0;
454 
455  while (it2.isNext()) {
456  auto p = it2.get();
457 
458  if (fabs(domain.getProp<1>(p) - domain.getProp<2>(p)) > worst) {
459  worst = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
460  }
461 
462  ++it2;
463  }
464  errv=v-vv;
465  domain.deleteGhost();
466  BOOST_REQUIRE(worst < 0.3);
467  }
468 
469  BOOST_AUTO_TEST_CASE(dcpse_op_div) {
470 // int rank;
471 // MPI_Comm_rank(MPI_COMM_WORLD, &rank);
472  size_t edgeSemiSize = 31;
473  const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
474  Box<2, double> box({0, 0}, {10.0, 10.0});
475  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
476  double spacing[2];
477  spacing[0] = box.getHigh(0) / (sz[0] - 1);
478  spacing[1] = box.getHigh(1) / (sz[1] - 1);
479  double rCut = 3.1* spacing[0];
480  Ghost<2, double> ghost(rCut);
481  BOOST_TEST_MESSAGE("Init vector_dist...");
482  double sigma2 = spacing[0] * spacing[1] / (2 * 4);
483 
484  vector_dist<2, double, aggregate<double, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>,double,double>> domain(
485  0, box, bc, ghost);
486 
487  //Init_DCPSE(domain)
488  BOOST_TEST_MESSAGE("Init domain...");
489 // std::random_device rd{};
490 // std::mt19937 rng{rd()};
491  std::mt19937 rng{6666666};
492 
493  std::normal_distribution<> gaussian{0, sigma2};
494 
495  auto it = domain.getGridIterator(sz);
496  size_t pointId = 0;
497  size_t counter = 0;
498  double minNormOne = 999;
499  while (it.isNext()) {
500  domain.add();
501  auto key = it.get();
502  mem_id k0 = key.get(0);
503  double x = k0 * spacing[0];
504  domain.getLastPos()[0] = x;//+ gaussian(rng);
505  mem_id k1 = key.get(1);
506  double y = k1 * spacing[1];
507  domain.getLastPos()[1] = y;//+gaussian(rng);
508  // Here fill the function value
509  domain.template getLastProp<1>()[0] = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
510  domain.template getLastProp<1>()[1] = cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]);
511 
512 
513  // Here fill the validation value for Divergence
514  domain.template getLastProp<0>()= cos(domain.getLastPos()[0]) - sin(domain.getLastPos()[1]);
515  /* domain.template getLastProp<4>()[0] =
516  cos(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) +
517  cos(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
518  domain.template getLastProp<4>()[1] =
519  -sin(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) -
520  sin(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));*/
521 
522 
523  ++counter;
524  ++it;
525  }
526  BOOST_TEST_MESSAGE("Sync domain across processors...");
527 
528  domain.map();
529  domain.ghost_get<0>();
530 
531  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
532  Divergence<decltype(verletList)> Div(domain, verletList, 2, rCut);
533  Derivative_x<decltype(verletList)> Dx(domain, verletList, 2, rCut);
534  Derivative_y<decltype(verletList)> Dy(domain, verletList, 2, rCut);
535 
536  auto v = getV<1>(domain);
537  auto anasol = getV<0>(domain);
538  auto div = getV<5>(domain);
539  auto div2 = getV<6>(domain);
540 
541  domain.ghost_get<1>();
542  div = Div(v);
543  div2=Dx(v[0])+Dy(v[1]);
544  auto it2 = domain.getDomainIterator();
545 
546  double worst1 = 0.0;
547  double worst2 = 0.0;
548 
549  while (it2.isNext()) {
550  auto p = it2.get();
551  if (fabs(domain.getProp<0>(p) - domain.getProp<5>(p)) > worst1) {
552  worst1 = fabs(domain.getProp<0>(p) - domain.getProp<5>(p));
553  }
554  if (fabs(domain.getProp<0>(p) - domain.getProp<6>(p)) > worst2) {
555  worst2 = fabs(domain.getProp<0>(p) - domain.getProp<6>(p));
556  }
557  ++it2;
558  }
559 
560  domain.deleteGhost();
561 /*
562  domain.write("DIV");
563 
564  std::cout<<worst1<<":"<<worst2;
565 */
566 
567  BOOST_REQUIRE(worst1 < 0.05);
568  BOOST_REQUIRE(worst2 < 0.05);
569 
570  }
571 
572  BOOST_AUTO_TEST_CASE(dcpse_op_vec) {
573 // int rank;
574 // MPI_Comm_rank(MPI_COMM_WORLD, &rank);
575  size_t edgeSemiSize = 160;
576  const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
577  Box<2, double> box({0, 0}, {10,10});
578  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
579  double spacing[2];
580  spacing[0] = box.getHigh(0) / (sz[0] - 1);
581  spacing[1] = box.getHigh(1) / (sz[1] - 1);
582  Ghost<2, double> ghost(spacing[0] * 3.1);
583  double rCut = 3.1 * spacing[0];
584  BOOST_TEST_MESSAGE("Init vector_dist...");
585  double sigma2 = spacing[0] * spacing[1] / (2 * 4);
586 
587  vector_dist<2, double, aggregate<double, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>,double,double>> domain(
588  0, box, bc, ghost);
589 
590  //Init_DCPSE(domain)
591  BOOST_TEST_MESSAGE("Init domain...");
592 // std::random_device rd{};
593 // std::mt19937 rng{rd()};
594  std::mt19937 rng{6666666};
595 
596  std::normal_distribution<> gaussian{0, sigma2};
597 
598  auto it = domain.getGridIterator(sz);
599  size_t pointId = 0;
600  size_t counter = 0;
601  double minNormOne = 999;
602  while (it.isNext()) {
603  domain.add();
604  auto key = it.get();
605  mem_id k0 = key.get(0);
606  double x = k0 * spacing[0];
607  domain.getLastPos()[0] = x;//+ gaussian(rng);
608  mem_id k1 = key.get(1);
609  double y = k1 * spacing[1];
610  domain.getLastPos()[1] = y;//+gaussian(rng);
611  // Here fill the function value
612  domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
613  domain.template getLastProp<1>()[0] = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
614  domain.template getLastProp<1>()[1] = cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]);
615 
616  // Here fill the validation value for Df/Dx
617  domain.template getLastProp<2>()[0] = 0;
618  domain.template getLastProp<2>()[1] = 0;
619  domain.template getLastProp<4>()[0] =
620  cos(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) +
621  cos(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
622  domain.template getLastProp<4>()[1] =
623  -sin(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) -
624  sin(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
625 
626  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]));
627 
628 
629  ++counter;
630  ++it;
631  }
632  BOOST_TEST_MESSAGE("Sync domain across processors...");
633 
634  domain.map();
635  domain.ghost_get<0>();
636 
637  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
638  Advection<decltype(verletList)> Adv(domain, verletList, 2, rCut);
639  auto v = getV<1>(domain);
640  auto P = getV<0>(domain);
641  auto dv = getV<3>(domain);
642  auto dP = getV<6>(domain);
643 
644 
645  domain.ghost_get<1>();
646  dv = Adv(v, v);
647  auto it2 = domain.getDomainIterator();
648 
649  double worst1 = 0.0;
650 
651  while (it2.isNext()) {
652  auto p = it2.get();
653 
654  if (fabs(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]) > worst1) {
655  worst1 = fabs(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]);
656 
657  }
658 
659  ++it2;
660  }
661 
662  //std::cout << "Maximum Error in component 2: " << worst1 << std::endl;
663  //domain.write("v1");
664  BOOST_REQUIRE(worst1 < 0.03);
665 
666 
667  dP = Adv(v, P);
668  auto it3 = domain.getDomainIterator();
669 
670  double worst2 = 0.0;
671 
672  while (it3.isNext()) {
673  auto p = it3.get();
674  if (fabs(domain.getProp<6>(p) - domain.getProp<5>(p)) > worst2) {
675  worst2 = fabs(domain.getProp<6>(p) - domain.getProp<5>(p));
676 
677  }
678 
679  ++it3;
680  }
681 
682 
683  domain.deleteGhost();
684  //domain.write("v2");
685  BOOST_REQUIRE(worst2 < 0.03);
686 
687  }
688 
689 
690  BOOST_AUTO_TEST_CASE(dcpse_slice) {
691  const size_t sz[2] = {31,31};
692  Box<2, double> box({0, 0}, {1,1});
693  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
694  double spacing = box.getHigh(0) / (sz[0] - 1);
695  Ghost<2, double> ghost(spacing * 3.9);
696  double rCut = 3.9 * spacing;
697 
698  vector_dist<2, double, aggregate<double,VectorS<2, double>,double,double[3][3]>> Particles(0, box, bc, ghost);
699 
700  auto it = Particles.getGridIterator(sz);
701  while (it.isNext()) {
702  Particles.add();
703  auto key = it.get();
704  double x = key.get(0) * it.getSpacing(0);
705  Particles.getLastPos()[0] = x;
706  double y = key.get(1) * it.getSpacing(1);
707  Particles.getLastPos()[1] = y;
708 
709  Particles.getLastProp<1>()[0] = sin(x+y);
710  Particles.getLastProp<1>()[1] = cos(x+y);
711 
712  ++it;
713  }
714 
715  Particles.map();
716  Particles.ghost_get<0,1>();
717 
718 
719  auto P = getV<0>(Particles);
720  auto V = getV<1>(Particles);
721  auto S = getV<2>(Particles);
722  auto Sig = getV<3>(Particles);
723 
724 
725  auto verletList = Particles.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
726  Derivative_x<decltype(verletList)> Dx(Particles, verletList, 2, rCut);
727 
728  P = Dx(V[0]);
729  S = V[0]*V[0] + V[1]*V[1];
730 
731  Sig[0][1] = V[0]*V[0] + V[1]*V[1];
732  Sig[1][0] = P;
733  Sig[2][2] = 5.0;
734 
735  auto it2 = Particles.getDomainIterator();
736 
737  double err = 0.0;
738 
739  while (it2.isNext())
740  {
741  auto p = it2.get();
742 
743  if (fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]) >= err )
744  {
745  err = fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]);
746  }
747 
748  if (fabs(Particles.getProp<2>(p) - 1.0) >= err )
749  {
750  err = fabs(Particles.getProp<2>(p) - 1.0);
751  }
752 
753  if (fabs(Particles.getProp<3>(p)[0][1] - 1.0) >= err )
754  {
755  err = fabs(Particles.getProp<3>(p)[0][1] - 1.0);
756  }
757 
758  if (fabs(Particles.getProp<3>(p)[1][0] - Particles.getProp<1>(p)[1]) >= err )
759  {
760  err = fabs(Particles.getProp<3>(p)[1][0] - Particles.getProp<1>(p)[1]);
761  }
762 
763  if (fabs(Particles.getProp<3>(p)[2][2] - 5.0) >= err )
764  {
765  err = fabs(Particles.getProp<3>(p)[2][2] - 5.0);
766  }
767 
768  ++it2;
769  }
770 
771  //Particles.write("test_out");
772  //std::cout << "Error: " << err << " " << create_vcluster().rank() << std::endl;
773  BOOST_REQUIRE(err < 0.03);
774 
775  }
776 
777  BOOST_AUTO_TEST_CASE(dcpse_slice_3d) {
778  const size_t sz[3] = {17,17,17};
779  Box<3, double> box({0, 0,0}, {1,1,1});
780  size_t bc[3] = {NON_PERIODIC, NON_PERIODIC,NON_PERIODIC};
781  double spacing = box.getHigh(0) / (sz[0] - 1);
782  Ghost<3, double> ghost(spacing * 3.9);
783  double rCut = 3.9 * spacing;
784 
785  vector_dist<3, double, aggregate<double,VectorS<3, double>,double,double[3][3]>> Particles(0, box, bc, ghost);
786 
787  auto it = Particles.getGridIterator(sz);
788  while (it.isNext()) {
789  Particles.add();
790  auto key = it.get();
791  double x = key.get(0) * it.getSpacing(0);
792  Particles.getLastPos()[0] = x;
793  double y = key.get(1) * it.getSpacing(1);
794  Particles.getLastPos()[1] = y;
795  double z = key.get(2) * it.getSpacing(2);
796  Particles.getLastPos()[2] = z;
797 
798  Particles.getLastProp<1>()[0] = sin(x+y);
799  Particles.getLastProp<1>()[1] = cos(x+y);
800  Particles.getLastProp<1>()[2] = 1.0;
801 
802  ++it;
803  }
804 
805  Particles.map();
806  Particles.ghost_get<0,1>();
807 
808 
809  auto P = getV<0>(Particles);
810  auto V = getV<1>(Particles);
811  auto S = getV<2>(Particles);
812  auto Sig = getV<3>(Particles);
813 
814  auto verletList = Particles.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
815  Derivative_x<decltype(verletList)> Dx(Particles, verletList, 2, rCut);
816 
817  P = Dx(V[0]);
818  S = V[0]*V[0] + V[1]*V[1]+V[2]*V[2];
819 
820  Sig[0][1] = V[0]*V[0] + V[1]*V[1]+V[2]*V[2];
821  Sig[1][0] = P;
822  Sig[2][2] = 5.0;
823 
824  auto it2 = Particles.getDomainIterator();
825 
826  double err1 = 0.0;
827  double err2 = 0.0;
828  double err3 = 0.0;
829  double err4 = 0.0;
830  double err5 = 0.0;
831 
832  while (it2.isNext())
833  {
834  auto p = it2.get();
835 
836  // 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)
837  if (fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]) >= err1 )
838  {
839  err1 = fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]);
840  }
841 
842  if (fabs(Particles.getProp<2>(p) - 2.0) >= err2 )
843  {
844  err2 = fabs(Particles.getProp<2>(p) - 2.0);
845  }
846 
847  if (fabs(Particles.getProp<3>(p)[0][1] - 2.0) >= err3 )
848  {
849  err3 = fabs(Particles.getProp<3>(p)[0][1] - 2.0);
850  }
851 
852  // V[0]*V[0] + V[1]*V[1]+V[2]*V[2] = 2 and
853  if (fabs(Particles.getProp<3>(p)[0][1] - Particles.getProp<2>(p)) >= err4 )
854  {
855  err4 = fabs(Particles.getProp<3>(p)[0][1] - Particles.getProp<2>(p));
856  }
857 
858  if (fabs(Particles.getProp<3>(p)[2][2] - 5.0) >= err5 )
859  {
860  err5 = fabs(Particles.getProp<3>(p)[2][2] - 5.0);
861  }
862 
863  ++it2;
864  }
865 
866  //std::cout << err1 << " " << err2 << " " << err3 << " " << err4 << " " << err5 << std::endl;
867  BOOST_REQUIRE(err1 < 0.08);
868  BOOST_REQUIRE(err2 < 0.03);
869  BOOST_REQUIRE(err3 < 0.03);
870  BOOST_REQUIRE(err4 < 0.03);
871  BOOST_REQUIRE(err5 < 0.03);
872  }
873 
874 
875 // Added by foggia on 08.03.2024
876 BOOST_AUTO_TEST_CASE(slicer_3index_tensor) {
877 
878  Vcluster<> & v_cl = create_vcluster();
879 
880  constexpr int VEC{0}; // vector - [27] linearized version of [DIM]x[DIM]x[DIM]
881  constexpr int MAT{1}; // matrix - [DIM]x[DIM]x[DIM]
882  constexpr int ERR{2};
883 
884  typedef aggregate<double[3*3*3],double[3][3][3],double[3*3*3]> prop_part;
885  openfpm::vector<std::string> propNames{"vector","matrix","err"};
886  typedef vector_dist<3,double,prop_part> vector_type;
887 
888  int n_part{1};
889  double rCut{3.0};
890  double grid_spacing{2.0/32.0};
891  Box<3,double> domain{{-1,-1,-1},{1,1,1}};
892  Ghost<3,double> ghost{rCut*grid_spacing + grid_spacing/8.0};
893  size_t bc[3] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
894 
895  vector_type part{0,domain,bc,ghost};
896  part.setPropNames(propNames);
897 
898  std::array<double,3> center{0,0,0};
899  std::array<double,3> coord;
900 
901  if (v_cl.rank() == 0) {
902  std::cout << "1. Creating particles\n";
903 
904  // Created using the Fibonacci sphere algorithm
905  const double pi{3.14159265358979323846};
906  const double golden_ang = pi*(3.0 - std::sqrt(5.0));
907  const double prefactor{std::sqrt(0.75/pi)};
908  double rad, theta, arg, thetaB, phi, phi_norm;
909 
910  for (int i = 0; i < n_part; ++i) {
911 
912  coord[1] = 1.0 - 2.0*(i/double(n_part-1));
913  rad = std::sqrt(1.0 - (coord[1]-center[1])*(coord[1]-center[1]));
914  theta = golden_ang * i;
915  coord[0] = std::cos(theta) * rad;
916  coord[2] = std::sin(theta) * rad;
917 
918  arg = (coord[0]-center[0]) * (coord[0]-center[0]) + (coord[1]-center[1]) * (coord[1]-center[1]);
919  thetaB = std::atan2(std::sqrt(arg),(coord[2]-center[2]));
920  phi = std::atan2((coord[1]-center[1]),(coord[0]-center[0]));
921 
922  part.add();
923  part.getLastPos()[0] = coord[0];
924  part.getLastPos()[1] = coord[1];
925  part.getLastPos()[2] = coord[2];
926 
927  for (int l = 0; l < 3; ++l)
928  for (int k = 0; k < 3; ++k)
929  for (int m = 0; m < 3; ++m) {
930 
931  int linIdx = (l*3 + k)*3 + m;
932 
933  part.getLastProp<VEC>()[linIdx] = double(linIdx);
934  part.getLastProp<MAT>()[l][k][m] = double(linIdx);
935  }
936  }
937  }
938 
939  part.map();
940  part.ghost_get<VEC>();
941 
942  double eps{1e-13};
943  auto vec{getV<VEC>(part)};
944  auto mat{getV<MAT>(part)};
945  auto err{getV<ERR>(part)};
946 
947  err[0] = mat[0][0][0] - vec[0];
948  err[1] = mat[0][0][1] - vec[1];
949  err[2] = mat[0][0][2] - vec[2];
950  err[3] = mat[0][1][0] - vec[3];
951  err[4] = mat[0][1][1] - vec[4];
952  err[5] = mat[0][1][2] - vec[5];
953  err[6] = mat[0][2][0] - vec[6];
954  err[7] = mat[0][2][1] - vec[7];
955  err[8] = mat[0][2][2] - vec[8];
956  err[9] = mat[1][0][0] - vec[9];
957  err[10] = mat[1][0][1] - vec[10];
958  err[11] = mat[1][0][2] - vec[11];
959  err[12] = mat[1][1][0] - vec[12];
960  err[13] = mat[1][1][1] - vec[13];
961  err[14] = mat[1][1][2] - vec[14];
962  err[15] = mat[1][2][0] - vec[15];
963  err[16] = mat[1][2][1] - vec[16];
964  err[17] = mat[1][2][2] - vec[17];
965  err[18] = mat[2][0][0] - vec[18];
966  err[19] = mat[2][0][1] - vec[19];
967  err[20] = mat[2][0][2] - vec[20];
968  err[21] = mat[2][1][0] - vec[21];
969  err[22] = mat[2][1][1] - vec[22];
970  err[23] = mat[2][1][2] - vec[23];
971  err[24] = mat[2][2][0] - vec[24];
972  err[25] = mat[2][2][1] - vec[25];
973  err[26] = mat[2][2][2] - vec[26];
974 
975  auto pit{part.getDomainIterator()};
976  while (pit.isNext()) {
977  auto key{pit.get()};
978 
979  // Uncomment to check the values
980  // for (int l = 0; l < 3; ++l)
981  // for (int k = 0; k < 3; ++k)
982  // for (int m = 0; m < 3; ++m)
983  // std::cout << "mat[" << l << "][" << k << "][" << m << "]: " << part.getProp<MAT>(key)[l][k][m] << std::endl;
984 
985  // for (int i = 0; i < 3*3*3; ++i) {
986  // std::cout << "vec[" << i << "]: " << part.getProp<VEC>(key)[i] << std::endl;
987  // std::cout << "err[" << i << "]: " << part.getProp<ERR>(key)[i] << std::endl;
988  // }
989 
990  for (int i = 0; i < 3*3*3; ++i) {
991  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[0],0, eps);
992  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[1],0, eps);
993  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[2],0, eps);
994  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[3],0, eps);
995  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[4],0, eps);
996  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[5],0, eps);
997  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[6],0, eps);
998  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[7],0, eps);
999  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[8],0, eps);
1000  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[9],0, eps);
1001  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[10],0, eps);
1002  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[11],0, eps);
1003  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[12],0, eps);
1004  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[13],0, eps);
1005  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[14],0, eps);
1006  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[15],0, eps);
1007  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[16],0, eps);
1008  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[17],0, eps);
1009  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[18],0, eps);
1010  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[19],0, eps);
1011  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[20],0, eps);
1012  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[21],0, eps);
1013  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[22],0, eps);
1014  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[23],0, eps);
1015  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[24],0, eps);
1016  BOOST_REQUIRE_CLOSE(part.getProp<ERR>(key)[25],0, eps);
1017  }
1018 
1019  ++pit;
1020  }
1021  part.deleteGhost();
1022  // part.write("test_3index");
1023 }
1024 
1025  BOOST_AUTO_TEST_CASE(dcpse_op_convection) {
1026  size_t edgeSemiSize = 20;
1027  const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
1028  Box<2, double> box({-1, -1}, {1.0, 1.0});
1029  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
1030  double spacing[2];
1031  spacing[0] = 2.0/ (sz[0] - 1);
1032  spacing[1] = 2.0 / (sz[1] - 1);
1033  double rCut = 3.1 * spacing[0];
1034  Ghost<2, double> ghost(rCut);
1035 
1036  BOOST_TEST_MESSAGE("Init vector_dist...");
1037 
1038  vector_dist<2, double, aggregate<double, double, double, VectorS<2, double>>> domain(0, box,bc,
1039  ghost);
1040  domain.setPropNames({"Concentration","Concentration_temp","Temp","Velocity"});
1041 
1042  //Init_DCPSE(domain)
1043  BOOST_TEST_MESSAGE("Init domain...");
1044 
1045  auto it = domain.getGridIterator(sz);
1046  while (it.isNext()) {
1047  domain.add();
1048  auto key = it.get();
1049  mem_id k0 = key.get(0);
1050  double x = -1.0+k0 * spacing[0];
1051  domain.getLastPos()[0] = x;//+ gaussian(rng);
1052  mem_id k1 = key.get(1);
1053  double y = -1.0+k1 * spacing[1];
1054  domain.getLastPos()[1] = y;//+gaussian(rng);
1055  // Here fill the function value
1056  if (x>-1 && y>-1 && x<1 && y<1)
1057  {
1058  domain.template getLastProp<3>()[0] = (-y)*exp(-10*((x)*(x)+(y)*(y)));;
1059  domain.template getLastProp<3>()[1] = (x)*exp(-10*((x)*(x)+(y)*(y)));;
1060  }
1061  else{
1062  domain.template getLastProp<3>()[0] = 0.0;
1063  domain.template getLastProp<3>()[1] = 0.0;
1064  }
1065  if (x==0.0 && y>-0.5 && y<0.5)
1066  {
1067  domain.template getLastProp<0>() = 1.0;
1068  }
1069  else
1070  {
1071  domain.template getLastProp<0>() = 0.0;
1072  }
1073  ++it;
1074  }
1075  BOOST_TEST_MESSAGE("Sync domain across processors...");
1076 
1077  domain.map();
1078  domain.ghost_get<0>();
1079 
1080  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
1081 
1082  //Derivative_x<decltype(verletList)> Dx(domain, verletList, 2, rCut);
1083  Derivative_xx<decltype(verletList)> Dxx(domain, verletList, 2, rCut);
1084  //Derivative_y<decltype(verletList)> Dy(domain, verletList, 2, rCut);
1085  Derivative_yy<decltype(verletList)> Dyy(domain, verletList, 2, rCut);
1086  auto C = getV<0>(domain);
1087  auto V = getV<3>(domain);
1088  auto Cnew = getV<1>(domain);
1089  auto Pos = getV<POS_PROP>(domain);
1090  timer tt;
1091  //domain.write_frame("Convection_init",0);
1092  int ctr=0;
1093  double t=0,tf=1,dt=1e-2;
1094  while(t<tf)
1095  {
1096  domain.write_frame("Convection",ctr);
1097  domain.ghost_get<0>();
1098  Cnew=C+dt*0.01*(Dxx(C)+Dyy(C));
1099  C=Cnew;
1100  Pos=Pos+dt*V;
1101  domain.map();
1102  domain.ghost_get<0>();
1103  auto it2 = domain.getDomainIterator();
1104  while (it2.isNext()) {
1105  auto p = it2.get();
1106  Point<2, double> xp = domain.getPos(p);
1107  double x=xp[0],y=xp[1];
1108  if (x>-1 && y>-1 && x<1 && y<1)
1109  {
1110  domain.template getProp<3>(p)[0] = (-y)*exp(-10*((x)*(x)+(y)*(y)));
1111  domain.template getProp<3>(p)[1] = (x)*exp(-10*((x)*(x)+(y)*(y)));;
1112  }
1113  else{
1114  domain.template getProp<3>(p)[0] = 0.0;
1115  domain.template getProp<3>(p)[1] = 0.0;
1116  }
1117  ++it2;
1118  }
1119  tt.start();
1120  domain.updateVerlet(verletList,rCut);
1121  Dxx.update(domain);
1122  Dyy.update(domain);
1123  tt.stop();
1124  //std::cout<tt.getwct()<<""
1125  ctr++;
1126  t+=dt;
1127  }
1128 
1129  Dxx.deallocate(domain);
1130  Dyy.deallocate(domain);
1131 
1132  /* std::cout<<"Dx"<<std::endl;
1133  Dx.checkMomenta(domain);
1134  std::cout<<"Dy"<<std::endl;
1135  Dy.checkMomenta(domain);
1136  std::cout<<"Dxx"<<std::endl;
1137  Dxx.checkMomenta(domain);
1138  int ctr=0;
1139  P=0;
1140  v=0;
1141  K1=0;
1142  auto its2 = domain.getDomainIterator();
1143  while (its2.isNext()) {
1144  auto p = its2.get();
1145  Dx.DrawKernel<0>(domain, p.getKey());
1146  Dy.DrawKernel<1>(domain, p.getKey());
1147  Dxx.DrawKernel<2>(domain, p.getKey());
1148  domain.write_frame("Kernels",ctr);
1149  P=0;
1150  v=0;
1151  K1=0;
1152  ++its2;
1153  ctr++;
1154  }*/
1155 
1156  }
1157 
1158 BOOST_AUTO_TEST_SUITE_END()
1159 
1160 
1161 #endif
1162 #endif
1163 
1164 
1165 
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
PPInterpolation
Class to perform particle to particle interpolation using DC-PSE kernels.
Definition: DcpseInterpolation.hpp:32
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
Vcluster_base::rank
size_t rank()
Get the process unit id.
Definition: VCluster_base.hpp:571
Vcluster
Implementation of VCluster class.
Definition: VCluster.hpp:59
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
aggregate
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Definition: aggregate.hpp:221