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