OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
DCPSE_op_test_base_tests.cu
1/*
2 * DCPSE_op_test_base_tests.cpp
3 *
4 * Created on: April 9, 2020
5 * Author: Abhinav Singh
6 *
7 */
8#include "config.h"
9#ifdef HAVE_EIGEN
10#ifdef HAVE_PETSC
11#define BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS
12#define BOOST_MPL_LIMIT_VECTOR_SIZE 30
13
14
15
16#define BOOST_TEST_DYN_LINK
17
18#include "util/util_debug.hpp"
19#include <boost/test/unit_test.hpp>
20#include <iostream>
21#include "../DCPSE_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_cu)
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_gpu<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>>> domain(0, box,
43 bc,
44 ghost);
45
46 //Init_DCPSE(domain)
47 BOOST_TEST_MESSAGE("Init domain...");
48
49 auto it = domain.getGridIterator(sz);
50 size_t pointId = 0;
51 size_t counter = 0;
52 double minNormOne = 999;
53 while (it.isNext()) {
54 domain.add();
55 auto key = it.get();
56 mem_id k0 = key.get(0);
57 double x = k0 * spacing[0];
58 domain.getLastPos()[0] = x;//+ gaussian(rng);
59 mem_id k1 = key.get(1);
60 double y = k1 * spacing[1];
61 domain.getLastPos()[1] = y;//+gaussian(rng);
62 // Here fill the function value
63 domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]);
64 domain.template getLastProp<2>() = 2*cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]);
65 ++counter;
66 ++it;
67 }
68 BOOST_TEST_MESSAGE("Sync domain across processors...");
69
70 domain.map();
71 domain.ghost_get<0>();
72 Derivative_x_gpu Dx(domain, 2, rCut);
73 Derivative_y_gpu Dy(domain, 2, rCut);
74 Gradient_gpu Grad(domain, 2, rCut);
75 Laplacian_gpu 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_gpu<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_gpu Dx(domain, 2, rCut);
143 Derivative_y_gpu 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_gpu DxLoaded(domain, 2, rCut,1,support_options::LOAD);
151 Derivative_y_gpu 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_gpu<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_gpu DxLoaded(domain, 2, rCut,1,support_options::LOAD);
218 Derivative_y_gpu 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_gpu<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_gpu<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>>> vector_type;
317
318 vector_type domain(0, box,bc,ghost);
319 vector_type 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_type,vector_type> 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_gpu<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_gpu 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_gpu<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_gpu Div(domain, 2, rCut);
519 Derivative_x_gpu Dx(domain, 2, rCut);
520 Derivative_y_gpu 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_gpu<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_gpu 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_gpu<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 auto P = getV<0>(Particles);
704 auto V = getV<1>(Particles);
705 auto S = getV<2>(Particles);
706 auto Sig = getV<3>(Particles);
707
708 Derivative_x_gpu Dx(Particles, 2, rCut,2);
709
710 P = Dx(V[0]);
711 S = V[0]*V[0] + V[1]*V[1];
712
713 Sig[0][1] = V[0]*V[0] + V[1]*V[1];
714 Sig[1][0] = P;
715 Sig[2][2] = 5.0;
716
717 auto it2 = Particles.getDomainIterator();
718
719 double err = 0.0;
720
721 while (it2.isNext())
722 {
723 auto p = it2.get();
724
725 if (fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]) >= err )
726 {
727 err = fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]);
728 }
729
730 if (fabs(Particles.getProp<2>(p) - 1.0) >= err )
731 {
732 err = fabs(Particles.getProp<2>(p) - 1.0);
733 }
734
735 if (fabs(Particles.getProp<3>(p)[0][1] - 1.0) >= err )
736 {
737 err = fabs(Particles.getProp<3>(p)[0][1] - 1.0);
738 }
739
740 if (fabs(Particles.getProp<3>(p)[1][0] - Particles.getProp<1>(p)[1]) >= err )
741 {
742 err = fabs(Particles.getProp<3>(p)[1][0] - Particles.getProp<1>(p)[1]);
743 }
744
745 if (fabs(Particles.getProp<3>(p)[2][2] - 5.0) >= err )
746 {
747 err = fabs(Particles.getProp<3>(p)[2][2] - 5.0);
748 }
749
750 ++it2;
751 }
752
753 //Particles.write("test_out");
754 //std::cout << "Error: " << err << " " << create_vcluster().rank() << std::endl;
755 BOOST_REQUIRE(err < 0.03);
756}
757
758 BOOST_AUTO_TEST_CASE(dcpse_slice_3d) {
759 const size_t sz[3] = {17,17,17};
760 Box<3, double> box({0, 0,0}, {1,1,1});
761 size_t bc[3] = {NON_PERIODIC, NON_PERIODIC,NON_PERIODIC};
762 double spacing = box.getHigh(0) / (sz[0] - 1);
763 Ghost<3, double> ghost(spacing * 3.9);
764 double rCut = 3.9 * spacing;
765
766 vector_dist_gpu<3, double, aggregate<double,VectorS<3, double>,double,double[3][3]>> Particles(0, box, bc, ghost);
767
768 auto it = Particles.getGridIterator(sz);
769 while (it.isNext()) {
770 Particles.add();
771 auto key = it.get();
772 double x = key.get(0) * it.getSpacing(0);
773 Particles.getLastPos()[0] = x;
774 double y = key.get(1) * it.getSpacing(1);
775 Particles.getLastPos()[1] = y;
776 double z = key.get(2) * it.getSpacing(2);
777 Particles.getLastPos()[2] = z;
778
779 Particles.getLastProp<1>()[0] = sin(x+y);
780 Particles.getLastProp<1>()[1] = cos(x+y);
781 Particles.getLastProp<1>()[2] = 1.0;
782
783 ++it;
784 }
785
786 Particles.map();
787 Particles.ghost_get<0,1>();
788
789
790 auto P = getV<0>(Particles);
791 auto V = getV<1>(Particles);
792 auto S = getV<2>(Particles);
793 auto Sig = getV<3>(Particles);
794
795
796 Derivative_x_gpu Dx(Particles, 2, rCut,2);
797
798 P = Dx(V[0]);
799 S = V[0]*V[0] + V[1]*V[1]+V[2]*V[2];
800
801 Sig[0][1] = V[0]*V[0] + V[1]*V[1]+V[2]*V[2];
802 Sig[1][0] = P;
803 Sig[2][2] = 5.0;
804
805 auto it2 = Particles.getDomainIterator();
806
807 double err1 = 0.0;
808 double err2 = 0.0;
809 double err3 = 0.0;
810 double err4 = 0.0;
811 double err5 = 0.0;
812
813 while (it2.isNext())
814 {
815 auto p = it2.get();
816
817 // 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)
818 if (fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]) >= err1 )
819 {
820 err1 = fabs(Particles.getProp<0>(p) - Particles.getProp<1>(p)[1]);
821 }
822
823 if (fabs(Particles.getProp<2>(p) - 2.0) >= err2 )
824 {
825 err2 = fabs(Particles.getProp<2>(p) - 2.0);
826 }
827
828 if (fabs(Particles.getProp<3>(p)[0][1] - 2.0) >= err3 )
829 {
830 err3 = fabs(Particles.getProp<3>(p)[0][1] - 2.0);
831 }
832
833 // V[0]*V[0] + V[1]*V[1]+V[2]*V[2] = 2 and
834 if (fabs(Particles.getProp<3>(p)[0][1] - Particles.getProp<2>(p)) >= err4 )
835 {
836 err4 = fabs(Particles.getProp<3>(p)[0][1] - Particles.getProp<2>(p));
837 }
838
839 if (fabs(Particles.getProp<3>(p)[2][2] - 5.0) >= err5 )
840 {
841 err5 = fabs(Particles.getProp<3>(p)[2][2] - 5.0);
842 }
843
844 ++it2;
845 }
846
847 //std::cout << err1 << " " << err2 << " " << err3 << " " << err4 << " " << err5 << std::endl;
848 BOOST_REQUIRE(err1 < 0.08);
849 BOOST_REQUIRE(err2 < 0.03);
850 BOOST_REQUIRE(err3 < 0.03);
851 BOOST_REQUIRE(err4 < 0.03);
852 BOOST_REQUIRE(err5 < 0.03);
853 }
854
855 BOOST_AUTO_TEST_CASE(dcpse_op_convection) {
856 size_t edgeSemiSize = 20;
857 const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
858 Box<2, double> box({-1, -1}, {1.0, 1.0});
859 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
860 double spacing[2];
861 spacing[0] = 2.0/ (sz[0] - 1);
862 spacing[1] = 2.0 / (sz[1] - 1);
863 double rCut = 3.1 * spacing[0];
864 Ghost<2, double> ghost(rCut);
865
866 BOOST_TEST_MESSAGE("Init vector_dist...");
867
868 vector_dist_gpu<2, double, aggregate<double, double, double, VectorS<2, double>>> domain(0, box,bc,
869 ghost);
870 domain.setPropNames({"Concentration","Concentration_temp","Temp","Velocity"});
871
872 //Init_DCPSE(domain)
873 BOOST_TEST_MESSAGE("Init domain...");
874
875 auto it = domain.getGridIterator(sz);
876 while (it.isNext()) {
877 domain.add();
878 auto key = it.get();
879 mem_id k0 = key.get(0);
880 double x = -1.0+k0 * spacing[0];
881 domain.getLastPos()[0] = x;//+ gaussian(rng);
882 mem_id k1 = key.get(1);
883 double y = -1.0+k1 * spacing[1];
884 domain.getLastPos()[1] = y;//+gaussian(rng);
885 // Here fill the function value
886 if (x>-1 && y>-1 && x<1 && y<1)
887 {
888 domain.template getLastProp<3>()[0] = (-y)*exp(-10*((x)*(x)+(y)*(y)));;
889 domain.template getLastProp<3>()[1] = (x)*exp(-10*((x)*(x)+(y)*(y)));;
890 }
891 else{
892 domain.template getLastProp<3>()[0] = 0.0;
893 domain.template getLastProp<3>()[1] = 0.0;
894 }
895 if (x==0.0 && y>-0.5 && y<0.5)
896 {
897 domain.template getLastProp<0>() = 1.0;
898 }
899 else
900 {
901 domain.template getLastProp<0>() = 0.0;
902 }
903 ++it;
904 }
905 BOOST_TEST_MESSAGE("Sync domain across processors...");
906
907 domain.map();
908 domain.ghost_get<0>();
909
910 //Derivative_x_gpu Dx(domain, 2, rCut);
911 Derivative_xx Dxx(domain, 2, rCut);
912 //Derivative_y_gpu Dy(domain, 2, rCut);
913 Derivative_yy Dyy(domain, 2, rCut);
914 auto C = getV<0>(domain);
915 auto V = getV<3>(domain);
916 auto Cnew = getV<1>(domain);
917 auto Pos = getV<PROP_POS>(domain);
918 timer tt;
919 //domain.write_frame("Convection_init",0);
920 int ctr=0;
921 double t=0,tf=1,dt=1e-2;
922 while(t<tf)
923 {
924 domain.write_frame("Convection",ctr);
925 domain.ghost_get<0>();
926 Cnew=C+dt*0.01*(Dxx(C)+Dyy(C));
927 C=Cnew;
928 Pos=Pos+dt*V;
929 domain.map();
930 domain.ghost_get<0>();
931 auto it2 = domain.getDomainIterator();
932 while (it2.isNext()) {
933 auto p = it2.get();
934 Point<2, double> xp = domain.getPos(p);
935 double x=xp[0],y=xp[1];
936 if (x>-1 && y>-1 && x<1 && y<1)
937 {
938 domain.template getProp<3>(p)[0] = (-y)*exp(-10*((x)*(x)+(y)*(y)));
939 domain.template getProp<3>(p)[1] = (x)*exp(-10*((x)*(x)+(y)*(y)));;
940 }
941 else{
942 domain.template getProp<3>(p)[0] = 0.0;
943 domain.template getProp<3>(p)[1] = 0.0;
944 }
945 ++it2;
946 }
947 tt.start();
948 Dxx.update(domain);
949 Dyy.update(domain);
950 tt.stop();
951 //std::cout<tt.getwct()<<""
952 ctr++;
953 t+=dt;
954 }
955
956 Dxx.deallocate(domain);
957 Dyy.deallocate(domain);
958
959 /* std::cout<<"Dx"<<std::endl;
960 Dx.checkMomenta(domain);
961 std::cout<<"Dy"<<std::endl;
962 Dy.checkMomenta(domain);
963 std::cout<<"Dxx"<<std::endl;
964 Dxx.checkMomenta(domain);
965 int ctr=0;
966 P=0;
967 v=0;
968 K1=0;
969 auto its2 = domain.getDomainIterator();
970 while (its2.isNext()) {
971 auto p = its2.get();
972 Dx.DrawKernel<0>(domain, p.getKey());
973 Dy.DrawKernel<1>(domain, p.getKey());
974 Dxx.DrawKernel<2>(domain, p.getKey());
975 domain.write_frame("Kernels",ctr);
976 P=0;
977 v=0;
978 K1=0;
979 ++its2;
980 ctr++;
981 }*/
982
983 }
984
985
986
987BOOST_AUTO_TEST_SUITE_END()
988
989
990#endif
991#endif
992
993
994
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