OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
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#include "DCPSE/DcpseInterpolation.hpp"
27
28BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests)
29BOOST_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<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 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 = 2*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 BOOST_AUTO_TEST_CASE(dcpse_op_save_load) {
100 size_t edgeSemiSize = 40;
101 const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
102 Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
103 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
104 double spacing[2];
105 spacing[0] = 2 * M_PI / (sz[0] - 1);
106 spacing[1] = 2 * M_PI / (sz[1] - 1);
107 Ghost<2, double> ghost(spacing[0] * 3.9);
108 double rCut = 3.9 * spacing[0];
109 BOOST_TEST_MESSAGE("Init vector_dist...");
110 double sigma2 = spacing[0] * spacing[1] / (2 * 4);
111
112 vector_dist<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>>> domain(0, box,
113 bc,
114 ghost);
115
116 //Init_DCPSE(domain)
117 BOOST_TEST_MESSAGE("Init domain...");
118
119 auto it = domain.getGridIterator(sz);
120 size_t pointId = 0;
121 size_t counter = 0;
122 double minNormOne = 999;
123 while (it.isNext()) {
124 domain.add();
125 auto key = it.get();
126 mem_id k0 = key.get(0);
127 double x = k0 * spacing[0];
128 domain.getLastPos()[0] = x;//+ gaussian(rng);
129 mem_id k1 = key.get(1);
130 double y = k1 * spacing[1];
131 domain.getLastPos()[1] = y;//+gaussian(rng);
132 // Here fill the function value
133 domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
134 domain.template getLastProp<2>() = 2*cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]);
135 ++counter;
136 ++it;
137 }
138 BOOST_TEST_MESSAGE("Sync domain across processors...");
139
140 domain.map();
141 domain.ghost_get<0>();
142 Derivative_x Dx(domain, 2, rCut);
143 Derivative_y Dy(domain, 2, rCut);
144 auto v = getV<1>(domain);
145 auto v2 = getV<3>(domain);
146 auto P = getV<0>(domain);
147 v2 = 2*Dx(P) + Dy(P);
148 Dx.save(domain,"DX_test");
149 Dy.save(domain,"DY_test");
150 Derivative_x DxLoaded(domain, 2, rCut,1,support_options::LOAD);
151 Derivative_y DyLoaded(domain, 2, rCut,1,support_options::LOAD);
152 DxLoaded.load(domain,"DX_test");
153 DyLoaded.load(domain,"DY_test");
154 v= 2*DxLoaded(P)+DyLoaded(P);
155 auto it2 = domain.getDomainIterator();
156 double worst = 0.0;
157 while (it2.isNext()) {
158 auto p = it2.get();
159
160 if (fabs(domain.getProp<1>(p) - domain.getProp<2>(p)) > worst) {
161 worst = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
162 }
163
164 ++it2;
165 }
166 domain.deleteGhost();
167 //std::cout<<worst;
168 BOOST_REQUIRE(worst < 0.03);
169 }
170
171 BOOST_AUTO_TEST_CASE(dcpse_op_save_load2) {
172 size_t edgeSemiSize = 40;
173 const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
174 Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
175 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
176 double spacing[2];
177 spacing[0] = 2 * M_PI / (sz[0] - 1);
178 spacing[1] = 2 * M_PI / (sz[1] - 1);
179 Ghost<2, double> ghost(spacing[0] * 3.9);
180 double rCut = 3.9 * spacing[0];
181 BOOST_TEST_MESSAGE("Init vector_dist...");
182 double sigma2 = spacing[0] * spacing[1] / (2 * 4);
183
184 vector_dist<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>>> domain(0, box,
185 bc,
186 ghost);
187
188 //Init_DCPSE(domain)
189 BOOST_TEST_MESSAGE("Init domain...");
190
191 auto it = domain.getGridIterator(sz);
192 size_t pointId = 0;
193 size_t counter = 0;
194 double minNormOne = 999;
195 while (it.isNext()) {
196 domain.add();
197 auto key = it.get();
198 mem_id k0 = key.get(0);
199 double x = k0 * spacing[0];
200 domain.getLastPos()[0] = x;//+ gaussian(rng);
201 mem_id k1 = key.get(1);
202 double y = k1 * spacing[1];
203 domain.getLastPos()[1] = y;//+gaussian(rng);
204 // Here fill the function value
205 domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
206 domain.template getLastProp<2>() = 2*cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]);
207 ++counter;
208 ++it;
209 }
210 BOOST_TEST_MESSAGE("Sync domain across processors...");
211
212 domain.map();
213 domain.ghost_get<0>();
214 auto v = getV<1>(domain);
215 auto v2 = getV<3>(domain);
216 auto P = getV<0>(domain);
217 Derivative_x DxLoaded(domain, 2, rCut,1,support_options::LOAD);
218 Derivative_y DyLoaded(domain, 2, rCut,1,support_options::LOAD);
219 DxLoaded.load(domain,"DX_test");
220 DyLoaded.load(domain,"DY_test");
221 v= 2*DxLoaded(P)+DyLoaded(P);
222 auto it2 = domain.getDomainIterator();
223 double worst = 0.0;
224 while (it2.isNext()) {
225 auto p = it2.get();
226
227 if (fabs(domain.getProp<1>(p) - domain.getProp<2>(p)) > worst) {
228 worst = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
229 }
230
231 ++it2;
232 }
233 domain.deleteGhost();
234 //std::cout<<worst;
235 BOOST_REQUIRE(worst < 0.03);
236 }
237
238 BOOST_AUTO_TEST_CASE(dcpse_op_tests_fa) {
239 size_t edgeSemiSize = 40;
240 const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
241 Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
242 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
243 double spacing[2];
244 spacing[0] = 2 * M_PI / (sz[0] - 1);
245 spacing[1] = 2 * M_PI / (sz[1] - 1);
246 Ghost<2, double> ghost(spacing[0] * 3.9);
247 double rCut = 3.9 * spacing[0];
248 BOOST_TEST_MESSAGE("Init vector_dist...");
249 double sigma2 = spacing[0] * spacing[1] / (2 * 4);
250
251 typedef vector_dist<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>>> vector_type;
252
253 vector_type domain(0, box,bc,ghost);
254
255 //Init_DCPSE(domain)
256 BOOST_TEST_MESSAGE("Init domain...");
257
258 auto it = domain.getGridIterator(sz);
259 size_t pointId = 0;
260 size_t counter = 0;
261 double minNormOne = 999;
262 while (it.isNext()) {
263 domain.add();
264 auto key = it.get();
265 mem_id k0 = key.get(0);
266 double x = k0 * spacing[0];
267 domain.getLastPos()[0] = x;//+ gaussian(rng);
268 mem_id k1 = key.get(1);
269 double y = k1 * spacing[1];
270 domain.getLastPos()[1] = y;//+gaussian(rng);
271 // Here fill the function value
272 domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
273 domain.template getLastProp<2>() = cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]);
274 ++counter;
275 ++it;
276 }
277 BOOST_TEST_MESSAGE("Sync domain across processors...");
278
279 domain.map();
280 domain.ghost_get<0>();
281
282 PPInterpolation<vector_type,vector_type> Fx(domain,domain, 2, rCut);
283 auto v = getV<1>(domain);
284 auto P = getV<0>(domain);
285
286 Fx.p2p<0,1>();
287 auto it2 = domain.getDomainIterator();
288 double worst = 0.0;
289 while (it2.isNext()) {
290 auto p = it2.get();
291 if (fabs(domain.getProp<1>(p) - domain.getProp<0>(p)) > worst) {
292 worst = fabs(domain.getProp<1>(p) - domain.getProp<0>(p));
293 }
294 ++it2;
295 }
296 //std::cout<<"Worst:"<<worst<<std::endl;
297 domain.deleteGhost();
298 //domain.write_frame("test",0,0.024,BINARY);
299 BOOST_REQUIRE(worst < 0.03);
300 }
301
302 BOOST_AUTO_TEST_CASE(dcpse_op_tests_mfa) {
303 size_t edgeSemiSize = 40;
304 const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
305 Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
306 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
307 double spacing[2];
308 spacing[0] = 2 * M_PI / (sz[0] - 1);
309 spacing[1] = 2 * M_PI / (sz[1] - 1);
310 Ghost<2, double> ghost(spacing[0] * 3.9);
311 double rCut = 3.9 * spacing[0];
312 BOOST_TEST_MESSAGE("Init vector_dist...");
313 double sigma2 = spacing[0] * spacing[1] / ( 4);
314 std::normal_distribution<> gaussian{0, sigma2};
315 std::mt19937 rng{6666666};
316 typedef vector_dist<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>>> vector_dist;
317
318 vector_dist domain(0, box,bc,ghost);
319 vector_dist domain2(domain.getDecomposition(),0);
320
321 //Init_DCPSE(domain)
322 BOOST_TEST_MESSAGE("Init domain...");
323
324 auto it = domain.getGridIterator(sz);
325 size_t pointId = 0;
326 size_t counter = 0;
327 double minNormOne = 999;
328 while (it.isNext()) {
329 domain.add();
330 domain2.add();
331 auto key = it.get();
332 mem_id k0 = key.get(0);
333 mem_id k1 = key.get(1);
334 double x = k0 * spacing[0];
335 double y = k1 * spacing[1];
336 domain.getLastPos()[0] = x;//+ gaussian(rng);
337 domain.getLastPos()[1] = y;//+gaussian(rng);
338 if(x!=0 && y!=0 && x!=box.getHigh(0) && y!=box.getHigh(1)){
339 domain2.getLastPos()[0] = x+ gaussian(rng);
340 domain2.getLastPos()[1] = y+ gaussian(rng);
341 }
342 else{
343 domain2.getLastPos()[0] = x;
344 domain2.getLastPos()[1] = y;
345 }
346 // Here fill the function value
347 domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
348 domain.template getLastProp<1>() = 0.0;
349 domain2.template getLastProp<0>() = sin(domain2.getLastPos()[0]) + sin(domain2.getLastPos()[1]);
350 ++counter;
351 ++it;
352 }
353 BOOST_TEST_MESSAGE("Sync domain across processors...");
354
355 domain.map();
356 domain2.map();
357 domain.ghost_get<0>();
358 domain2.ghost_get<0>();
359
360 PPInterpolation<vector_dist,vector_dist> Fx(domain2,domain, 2, rCut);
361 //auto v = getV<1>(domain);
362 //auto P = getV<0>(domain);
363 Fx.p2p<0,1>();
364 auto it2 = domain.getDomainIterator();
365 double worst = 0.0;
366 while (it2.isNext()) {
367 auto p = it2.get();
368 //domain.template getProp<2>(p) = domain.getProp<1>(p) - domain.getProp<0>(p);
369 if (fabs(domain.getProp<1>(p) - domain.getProp<0>(p)) > worst) {
370 worst = fabs(domain.getProp<1>(p) - domain.getProp<0>(p));
371 }
372 ++it2;
373 }
374 //std::cout<<"Worst:"<<worst<<std::endl;
375 domain.deleteGhost();
376 //domain.write("test1");
377 //domain2.write("test2");
378 BOOST_REQUIRE(worst < 0.03);
379 }
380
381
382 BOOST_AUTO_TEST_CASE(dcpse_op_test_lap) {
383 size_t edgeSemiSize = 81;
384 const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
385 Box<2, double> box({0, 0}, {2 * M_PI, 2 * M_PI});
386 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
387 double spacing[2];
388 spacing[0] = 2 * M_PI / (sz[0] - 1);
389 spacing[1] = 2 * M_PI / (sz[1] - 1);
390 Ghost<2, double> ghost(spacing[0] * 3.9);
391 double rCut = 3.9 * spacing[0];
392 BOOST_TEST_MESSAGE("Init vector_dist...");
393 double sigma2 = spacing[0] * spacing[1] / (2 * 4);
394
395 vector_dist<2, double, aggregate<double, double, double, double, VectorS<2, double>>> domain(0, box,
396 bc,
397 ghost);
398
399 //Init_DCPSE(domain)
400 BOOST_TEST_MESSAGE("Init domain...");
401 std::normal_distribution<> gaussian{0, sigma2};
402
403 auto it = domain.getGridIterator(sz);
404 size_t pointId = 0;
405 size_t counter = 0;
406 double minNormOne = 999;
407 while (it.isNext()) {
408 domain.add();
409 auto key = it.get();
410 mem_id k0 = key.get(0);
411 double x = k0 * spacing[0];
412 domain.getLastPos()[0] = x;//+ gaussian(rng);
413 mem_id k1 = key.get(1);
414 double y = k1 * spacing[1];
415 domain.getLastPos()[1] = y;//+gaussian(rng);
416 // Here fill the function value
417 domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
418 // Here fill the validation value for Lap
419 domain.template getLastProp<1>() = - sin(domain.getLastPos()[0]) - sin(domain.getLastPos()[1]);
420 domain.template getLastProp<4>()[0] = -sin(domain.getLastPos()[0]);
421 domain.template getLastProp<4>()[1] = -sin(domain.getLastPos()[1]);
422
423 ++counter;
424 ++it;
425 }
426 BOOST_TEST_MESSAGE("Sync domain across processors...");
427
428 domain.map();
429 domain.ghost_get<0>();
430
431 Laplacian Lap(domain, 2, rCut);
432 auto v = getV<1>(domain);
433 auto P = getV<0>(domain);
434 auto vv = getV<2>(domain);
435 auto errv = getV<3>(domain);
436
437 vv = Lap(P);
438 auto it2 = domain.getDomainIterator();
439
440 double worst = 0.0;
441
442 while (it2.isNext()) {
443 auto p = it2.get();
444
445 if (fabs(domain.getProp<1>(p) - domain.getProp<2>(p)) > worst) {
446 worst = fabs(domain.getProp<1>(p) - domain.getProp<2>(p));
447 }
448
449 ++it2;
450 }
451 errv=v-vv;
452 domain.deleteGhost();
453 BOOST_REQUIRE(worst < 0.3);
454 }
455
456 BOOST_AUTO_TEST_CASE(dcpse_op_div) {
457// int rank;
458// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
459 size_t edgeSemiSize = 31;
460 const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
461 Box<2, double> box({0, 0}, {10.0, 10.0});
462 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
463 double spacing[2];
464 spacing[0] = box.getHigh(0) / (sz[0] - 1);
465 spacing[1] = box.getHigh(1) / (sz[1] - 1);
466 double rCut = 3.1* spacing[0];
467 Ghost<2, double> ghost(rCut);
468 BOOST_TEST_MESSAGE("Init vector_dist...");
469 double sigma2 = spacing[0] * spacing[1] / (2 * 4);
470
471 vector_dist<2, double, aggregate<double, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>,double,double>> domain(
472 0, box, bc, ghost);
473
474 //Init_DCPSE(domain)
475 BOOST_TEST_MESSAGE("Init domain...");
476// std::random_device rd{};
477// std::mt19937 rng{rd()};
478 std::mt19937 rng{6666666};
479
480 std::normal_distribution<> gaussian{0, sigma2};
481
482 auto it = domain.getGridIterator(sz);
483 size_t pointId = 0;
484 size_t counter = 0;
485 double minNormOne = 999;
486 while (it.isNext()) {
487 domain.add();
488 auto key = it.get();
489 mem_id k0 = key.get(0);
490 double x = k0 * spacing[0];
491 domain.getLastPos()[0] = x;//+ gaussian(rng);
492 mem_id k1 = key.get(1);
493 double y = k1 * spacing[1];
494 domain.getLastPos()[1] = y;//+gaussian(rng);
495 // Here fill the function value
496 domain.template getLastProp<1>()[0] = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
497 domain.template getLastProp<1>()[1] = cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]);
498
499
500 // Here fill the validation value for Divergence
501 domain.template getLastProp<0>()= cos(domain.getLastPos()[0]) - sin(domain.getLastPos()[1]);
502 /* domain.template getLastProp<4>()[0] =
503 cos(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) +
504 cos(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
505 domain.template getLastProp<4>()[1] =
506 -sin(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) -
507 sin(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));*/
508
509
510 ++counter;
511 ++it;
512 }
513 BOOST_TEST_MESSAGE("Sync domain across processors...");
514
515 domain.map();
516 domain.ghost_get<0>();
517
518 Divergence Div(domain, 2, rCut);
519 Derivative_x Dx(domain, 2, rCut);
520 Derivative_y Dy(domain, 2, rCut);
521
522 auto v = getV<1>(domain);
523 auto anasol = getV<0>(domain);
524 auto div = getV<5>(domain);
525 auto div2 = getV<6>(domain);
526
527 domain.ghost_get<1>();
528 div = Div(v);
529 div2=Dx(v[0])+Dy(v[1]);
530 auto it2 = domain.getDomainIterator();
531
532 double worst1 = 0.0;
533 double worst2 = 0.0;
534
535 while (it2.isNext()) {
536 auto p = it2.get();
537 if (fabs(domain.getProp<0>(p) - domain.getProp<5>(p)) > worst1) {
538 worst1 = fabs(domain.getProp<0>(p) - domain.getProp<5>(p));
539 }
540 if (fabs(domain.getProp<0>(p) - domain.getProp<6>(p)) > worst2) {
541 worst2 = fabs(domain.getProp<0>(p) - domain.getProp<6>(p));
542 }
543 ++it2;
544 }
545
546 domain.deleteGhost();
547/*
548 domain.write("DIV");
549
550 std::cout<<worst1<<":"<<worst2;
551*/
552
553 BOOST_REQUIRE(worst1 < 0.05);
554 BOOST_REQUIRE(worst2 < 0.05);
555
556 }
557
558 BOOST_AUTO_TEST_CASE(dcpse_op_vec) {
559// int rank;
560// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
561 size_t edgeSemiSize = 160;
562 const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
563 Box<2, double> box({0, 0}, {10,10});
564 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
565 double spacing[2];
566 spacing[0] = box.getHigh(0) / (sz[0] - 1);
567 spacing[1] = box.getHigh(1) / (sz[1] - 1);
568 Ghost<2, double> ghost(spacing[0] * 3.1);
569 double rCut = 3.1 * spacing[0];
570 BOOST_TEST_MESSAGE("Init vector_dist...");
571 double sigma2 = spacing[0] * spacing[1] / (2 * 4);
572
573 vector_dist<2, double, aggregate<double, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>, VectorS<2, double>,double,double>> domain(
574 0, box, bc, ghost);
575
576 //Init_DCPSE(domain)
577 BOOST_TEST_MESSAGE("Init domain...");
578// std::random_device rd{};
579// std::mt19937 rng{rd()};
580 std::mt19937 rng{6666666};
581
582 std::normal_distribution<> gaussian{0, sigma2};
583
584 auto it = domain.getGridIterator(sz);
585 size_t pointId = 0;
586 size_t counter = 0;
587 double minNormOne = 999;
588 while (it.isNext()) {
589 domain.add();
590 auto key = it.get();
591 mem_id k0 = key.get(0);
592 double x = k0 * spacing[0];
593 domain.getLastPos()[0] = x;//+ gaussian(rng);
594 mem_id k1 = key.get(1);
595 double y = k1 * spacing[1];
596 domain.getLastPos()[1] = y;//+gaussian(rng);
597 // Here fill the function value
598 domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
599 domain.template getLastProp<1>()[0] = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
600 domain.template getLastProp<1>()[1] = cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]);
601
602 // Here fill the validation value for Df/Dx
603 domain.template getLastProp<2>()[0] = 0;
604 domain.template getLastProp<2>()[1] = 0;
605 domain.template getLastProp<4>()[0] =
606 cos(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) +
607 cos(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
608 domain.template getLastProp<4>()[1] =
609 -sin(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) -
610 sin(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
611
612 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]));
613
614
615 ++counter;
616 ++it;
617 }
618 BOOST_TEST_MESSAGE("Sync domain across processors...");
619
620 domain.map();
621 domain.ghost_get<0>();
622
623 Advection Adv(domain, 2, rCut);
624 auto v = getV<1>(domain);
625 auto P = getV<0>(domain);
626 auto dv = getV<3>(domain);
627 auto dP = getV<6>(domain);
628
629
630 domain.ghost_get<1>();
631 dv = Adv(v, v);
632 auto it2 = domain.getDomainIterator();
633
634 double worst1 = 0.0;
635
636 while (it2.isNext()) {
637 auto p = it2.get();
638
639 if (fabs(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]) > worst1) {
640 worst1 = fabs(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]);
641
642 }
643
644 ++it2;
645 }
646
647 //std::cout << "Maximum Error in component 2: " << worst1 << std::endl;
648 //domain.write("v1");
649 BOOST_REQUIRE(worst1 < 0.03);
650
651
652 dP = Adv(v, P);
653 auto it3 = domain.getDomainIterator();
654
655 double worst2 = 0.0;
656
657 while (it3.isNext()) {
658 auto p = it3.get();
659 if (fabs(domain.getProp<6>(p) - domain.getProp<5>(p)) > worst2) {
660 worst2 = fabs(domain.getProp<6>(p) - domain.getProp<5>(p));
661
662 }
663
664 ++it3;
665 }
666
667
668 domain.deleteGhost();
669 //domain.write("v2");
670 BOOST_REQUIRE(worst2 < 0.03);
671
672 }
673
674
675 BOOST_AUTO_TEST_CASE(dcpse_slice) {
676 const size_t sz[2] = {31,31};
677 Box<2, double> box({0, 0}, {1,1});
678 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
679 double spacing = box.getHigh(0) / (sz[0] - 1);
680 Ghost<2, double> ghost(spacing * 3.9);
681 double rCut = 3.9 * spacing;
682
683 vector_dist<2, double, aggregate<double,VectorS<2, double>,double,double[3][3]>> Particles(0, box, bc, ghost);
684
685 auto it = Particles.getGridIterator(sz);
686 while (it.isNext()) {
687 Particles.add();
688 auto key = it.get();
689 double x = key.get(0) * it.getSpacing(0);
690 Particles.getLastPos()[0] = x;
691 double y = key.get(1) * it.getSpacing(1);
692 Particles.getLastPos()[1] = y;
693
694 Particles.getLastProp<1>()[0] = sin(x+y);
695 Particles.getLastProp<1>()[1] = cos(x+y);
696
697 ++it;
698 }
699
700 Particles.map();
701 Particles.ghost_get<0,1>();
702
703
704 auto P = getV<0>(Particles);
705 auto V = getV<1>(Particles);
706 auto S = getV<2>(Particles);
707 auto Sig = getV<3>(Particles);
708
709
710 Derivative_x Dx(Particles, 2, rCut,2);
711
712 P = Dx(V[0]);
713 S = V[0]*V[0] + V[1]*V[1];
714
715 Sig[0][1] = V[0]*V[0] + V[1]*V[1];
716 Sig[1][0] = P;
717 Sig[2][2] = 5.0;
718
719 auto it2 = Particles.getDomainIterator();
720
721 double err = 0.0;
722
723 while (it2.isNext())
724 {
725 auto p = it2.get();
726
727 if (fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]) >= err )
728 {
729 err = fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]);
730 }
731
732 if (fabs(Particles.getProp<2>(p) - 1.0) >= err )
733 {
734 err = fabs(Particles.getProp<2>(p) - 1.0);
735 }
736
737 if (fabs(Particles.getProp<3>(p)[0][1] - 1.0) >= err )
738 {
739 err = fabs(Particles.getProp<3>(p)[0][1] - 1.0);
740 }
741
742 if (fabs(Particles.getProp<3>(p)[1][0] - Particles.getProp<1>(p)[1]) >= err )
743 {
744 err = fabs(Particles.getProp<3>(p)[1][0] - Particles.getProp<1>(p)[1]);
745 }
746
747 if (fabs(Particles.getProp<3>(p)[2][2] - 5.0) >= err )
748 {
749 err = fabs(Particles.getProp<3>(p)[2][2] - 5.0);
750 }
751
752 ++it2;
753 }
754
755 //Particles.write("test_out");
756 //std::cout << "Error: " << err << " " << create_vcluster().rank() << std::endl;
757 BOOST_REQUIRE(err < 0.03);
758
759 }
760
761 BOOST_AUTO_TEST_CASE(dcpse_slice_3d) {
762 const size_t sz[3] = {17,17,17};
763 Box<3, double> box({0, 0,0}, {1,1,1});
764 size_t bc[3] = {NON_PERIODIC, NON_PERIODIC,NON_PERIODIC};
765 double spacing = box.getHigh(0) / (sz[0] - 1);
766 Ghost<3, double> ghost(spacing * 3.9);
767 double rCut = 3.9 * spacing;
768
769 vector_dist<3, double, aggregate<double,VectorS<3, double>,double,double[3][3]>> Particles(0, box, bc, ghost);
770
771 auto it = Particles.getGridIterator(sz);
772 while (it.isNext()) {
773 Particles.add();
774 auto key = it.get();
775 double x = key.get(0) * it.getSpacing(0);
776 Particles.getLastPos()[0] = x;
777 double y = key.get(1) * it.getSpacing(1);
778 Particles.getLastPos()[1] = y;
779 double z = key.get(2) * it.getSpacing(2);
780 Particles.getLastPos()[2] = z;
781
782 Particles.getLastProp<1>()[0] = sin(x+y);
783 Particles.getLastProp<1>()[1] = cos(x+y);
784 Particles.getLastProp<1>()[2] = 1.0;
785
786 ++it;
787 }
788
789 Particles.map();
790 Particles.ghost_get<0,1>();
791
792
793 auto P = getV<0>(Particles);
794 auto V = getV<1>(Particles);
795 auto S = getV<2>(Particles);
796 auto Sig = getV<3>(Particles);
797
798
799 Derivative_x Dx(Particles, 2, rCut,2);
800
801 P = Dx(V[0]);
802 S = V[0]*V[0] + V[1]*V[1]+V[2]*V[2];
803
804 Sig[0][1] = V[0]*V[0] + V[1]*V[1]+V[2]*V[2];
805 Sig[1][0] = P;
806 Sig[2][2] = 5.0;
807
808 auto it2 = Particles.getDomainIterator();
809
810 double err1 = 0.0;
811 double err2 = 0.0;
812 double err3 = 0.0;
813 double err4 = 0.0;
814 double err5 = 0.0;
815
816 while (it2.isNext())
817 {
818 auto p = it2.get();
819
820 // 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)
821 if (fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]) >= err1 )
822 {
823 err1 = fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]);
824 }
825
826 if (fabs(Particles.getProp<2>(p) - 2.0) >= err2 )
827 {
828 err2 = fabs(Particles.getProp<2>(p) - 2.0);
829 }
830
831 if (fabs(Particles.getProp<3>(p)[0][1] - 2.0) >= err3 )
832 {
833 err3 = fabs(Particles.getProp<3>(p)[0][1] - 2.0);
834 }
835
836 // V[0]*V[0] + V[1]*V[1]+V[2]*V[2] = 2 and
837 if (fabs(Particles.getProp<3>(p)[0][1] - Particles.getProp<2>(p)) >= err4 )
838 {
839 err4 = fabs(Particles.getProp<3>(p)[0][1] - Particles.getProp<2>(p));
840 }
841
842 if (fabs(Particles.getProp<3>(p)[2][2] - 5.0) >= err5 )
843 {
844 err5 = fabs(Particles.getProp<3>(p)[2][2] - 5.0);
845 }
846
847 ++it2;
848 }
849
850 //std::cout << err1 << " " << err2 << " " << err3 << " " << err4 << " " << err5 << std::endl;
851 BOOST_REQUIRE(err1 < 0.08);
852 BOOST_REQUIRE(err2 < 0.03);
853 BOOST_REQUIRE(err3 < 0.03);
854 BOOST_REQUIRE(err4 < 0.03);
855 BOOST_REQUIRE(err5 < 0.03);
856 }
857
858 BOOST_AUTO_TEST_CASE(dcpse_op_convection) {
859 size_t edgeSemiSize = 20;
860 const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
861 Box<2, double> box({-1, -1}, {1.0, 1.0});
862 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
863 double spacing[2];
864 spacing[0] = 2.0/ (sz[0] - 1);
865 spacing[1] = 2.0 / (sz[1] - 1);
866 double rCut = 3.1 * spacing[0];
867 Ghost<2, double> ghost(rCut);
868
869 BOOST_TEST_MESSAGE("Init vector_dist...");
870
871 vector_dist<2, double, aggregate<double, double, double, VectorS<2, double>>> domain(0, box,bc,
872 ghost);
873 domain.setPropNames({"Concentration","Concentration_temp","Temp","Velocity"});
874
875 //Init_DCPSE(domain)
876 BOOST_TEST_MESSAGE("Init domain...");
877
878 auto it = domain.getGridIterator(sz);
879 while (it.isNext()) {
880 domain.add();
881 auto key = it.get();
882 mem_id k0 = key.get(0);
883 double x = -1.0+k0 * spacing[0];
884 domain.getLastPos()[0] = x;//+ gaussian(rng);
885 mem_id k1 = key.get(1);
886 double y = -1.0+k1 * spacing[1];
887 domain.getLastPos()[1] = y;//+gaussian(rng);
888 // Here fill the function value
889 if (x>-1 && y>-1 && x<1 && y<1)
890 {
891 domain.template getLastProp<3>()[0] = (-y)*exp(-10*((x)*(x)+(y)*(y)));;
892 domain.template getLastProp<3>()[1] = (x)*exp(-10*((x)*(x)+(y)*(y)));;
893 }
894 else{
895 domain.template getLastProp<3>()[0] = 0.0;
896 domain.template getLastProp<3>()[1] = 0.0;
897 }
898 if (x==0.0 && y>-0.5 && y<0.5)
899 {
900 domain.template getLastProp<0>() = 1.0;
901 }
902 else
903 {
904 domain.template getLastProp<0>() = 0.0;
905 }
906 ++it;
907 }
908 BOOST_TEST_MESSAGE("Sync domain across processors...");
909
910 domain.map();
911 domain.ghost_get<0>();
912
913 //Derivative_x Dx(domain, 2, rCut);
914 Derivative_xx Dxx(domain, 2, rCut);
915 //Derivative_y Dy(domain, 2, rCut);
916 Derivative_yy Dyy(domain, 2, rCut);
917 auto C = getV<0>(domain);
918 auto V = getV<3>(domain);
919 auto Cnew = getV<1>(domain);
920 auto Pos = getV<PROP_POS>(domain);
921 timer tt;
922 //domain.write_frame("Convection_init",0);
923 int ctr=0;
924 double t=0,tf=1,dt=1e-2;
925 while(t<tf)
926 {
927 domain.write_frame("Convection",ctr);
928 domain.ghost_get<0>();
929 Cnew=C+dt*0.01*(Dxx(C)+Dyy(C));
930 C=Cnew;
931 Pos=Pos+dt*V;
932 domain.map();
933 domain.ghost_get<0>();
934 auto it2 = domain.getDomainIterator();
935 while (it2.isNext()) {
936 auto p = it2.get();
937 Point<2, double> xp = domain.getPos(p);
938 double x=xp[0],y=xp[1];
939 if (x>-1 && y>-1 && x<1 && y<1)
940 {
941 domain.template getProp<3>(p)[0] = (-y)*exp(-10*((x)*(x)+(y)*(y)));
942 domain.template getProp<3>(p)[1] = (x)*exp(-10*((x)*(x)+(y)*(y)));;
943 }
944 else{
945 domain.template getProp<3>(p)[0] = 0.0;
946 domain.template getProp<3>(p)[1] = 0.0;
947 }
948 ++it2;
949 }
950 tt.start();
951 Dxx.update(domain);
952 Dyy.update(domain);
953 tt.stop();
954 //std::cout<tt.getwct()<<""
955 ctr++;
956 t+=dt;
957 }
958
959 Dxx.deallocate(domain);
960 Dyy.deallocate(domain);
961
962 /* std::cout<<"Dx"<<std::endl;
963 Dx.checkMomenta(domain);
964 std::cout<<"Dy"<<std::endl;
965 Dy.checkMomenta(domain);
966 std::cout<<"Dxx"<<std::endl;
967 Dxx.checkMomenta(domain);
968 int ctr=0;
969 P=0;
970 v=0;
971 K1=0;
972 auto its2 = domain.getDomainIterator();
973 while (its2.isNext()) {
974 auto p = its2.get();
975 Dx.DrawKernel<0>(domain, p.getKey());
976 Dy.DrawKernel<1>(domain, p.getKey());
977 Dxx.DrawKernel<2>(domain, p.getKey());
978 domain.write_frame("Kernels",ctr);
979 P=0;
980 v=0;
981 K1=0;
982 ++its2;
983 ctr++;
984 }*/
985
986 }
987
988
989
990BOOST_AUTO_TEST_SUITE_END()
991
992
993#endif
994#endif
995
996
997
Box
This class represent an N-dimensional box.
Definition Box.hpp:61
Lap
Laplacian second order on h (spacing)
Definition Laplacian.hpp:23
PPInterpolation
Class for Creating the DCPSE Operator For the function approximation objects and computes DCPSE Kerne...
Definition DcpseInterpolation.hpp:23
Point_test
Test structure used for several test.
Definition Point_test.hpp:106
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
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:305