OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
DCPSE_op_Solver_test.cu
1/*
2 * DCPSE_op_Solver_test.cu
3 *
4 * Created on: Jan 7, 2020
5 * Author: Abhinav Singh, Pietro Incardona, Serhii
6 *
7 */
8#define BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS
9#define BOOST_MPL_LIMIT_VECTOR_SIZE 40
10#include "config.h"
11
12
13#define BOOST_TEST_DYN_LINK
14
15#include "util/util_debug.hpp"
16#include <boost/test/unit_test.hpp>
17#include <iostream>
18#include "../DCPSE_op.hpp"
19#include "../DCPSE_Solver.hpp"
20#include "Operators/Vector/vector_dist_operators.hpp"
21#include "Vector/vector_dist_subset.hpp"
22#include "../EqnsStruct.hpp"
23#include "Decomposition/Distribution/SpaceDistribution.hpp"
24
25BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests_cu)
26
27 BOOST_AUTO_TEST_CASE(dcpse_op_solver) {
28// int rank;
29// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
30 const size_t sz[2] = {31, 31};
31 Box<2, double> box({0, 0}, {1.0, 1.0});
32 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
33 double spacing = box.getHigh(0) / (sz[0] - 1);
34 Ghost<2, double> ghost(spacing * 3);
35 double rCut = 2.0 * spacing;
36 BOOST_TEST_MESSAGE("Init vector_dist...");
37
38 vector_dist_gpu<2, double, aggregate<double,double,double,double>> domain(0, box, bc, ghost);
39
40
41 //Init_DCPSE(domain)
42 BOOST_TEST_MESSAGE("Init domain...");
43
44 auto it = domain.getGridIterator(sz);
45 while (it.isNext()) {
46 domain.add();
47
48 auto key = it.get();
49 double x = key.get(0) * it.getSpacing(0);
50 domain.getLastPos()[0] = x;
51 double y = key.get(1) * it.getSpacing(1);
52 domain.getLastPos()[1] = y;
53
54 ++it;
55 }
56 BOOST_TEST_MESSAGE("Sync domain across processors...");
57
58 domain.map();
59 domain.ghost_get<0>();
60
61 Laplacian_gpu Lap(domain, 2, rCut, 2,support_options::N_PARTICLES);
62
63 DCPSE_scheme_gpu<equations2d1_gpu,decltype(domain)> Solver( domain);
64
65 openfpm::vector<aggregate<int>> bulk;
66 openfpm::vector<aggregate<int>> up_p;
67 openfpm::vector<aggregate<int>> dw_p;
68 openfpm::vector<aggregate<int>> l_p;
69 openfpm::vector<aggregate<int>> r_p;
70
71 auto v = getV<0>(domain);
72 auto RHS = getV<1>(domain);
73 auto sol = getV<2>(domain);
74 auto anasol = getV<3>(domain);
75
76 // Here fill me
77
78 Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
79 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
80
81 Box<2, double> up_d({box.getLow(0) - spacing / 2.0, box.getHigh(1) - 8*spacing / 2.0},
82 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - 6*spacing / 2.0});
83
84 Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
85 {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
86
87 Box<2, double> down_u({box.getLow(0) - spacing / 2.0, box.getLow(1) + 3*spacing / 2.0},
88 {box.getHigh(0) + spacing / 2.0, box.getLow(1) + 4*spacing / 2.0});
89
90 Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
91 {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
92
93 Box<2, double> left_r({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
94 {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
95
96 Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
97 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
98
99 Box<2, double> right_l({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
100 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
101
102
103 openfpm::vector<Box<2, double>> boxes;
104 boxes.add(up);
105 boxes.add(up_d);
106 boxes.add(down);
107 boxes.add(down_u);
108 boxes.add(left);
109 boxes.add(left_r);
110 boxes.add(right);
111 boxes.add(right_l);
112
113 // Create a writer and write
114 VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
115 vtk_box.add(boxes);
116 vtk_box.write("vtk_box.vtk");
117
118 auto it2 = domain.getDomainIterator();
119
120 while (it2.isNext()) {
121 auto p = it2.get();
122 Point<2, double> xp = domain.getPos(p);
123 domain.getProp<2>(p)=1+xp[0]*xp[0]+2*xp[1]*xp[1];
124 if (up.isInside(xp) == true) {
125 up_p.add();
126 up_p.last().get<0>() = p.getKey();
127 domain.getProp<1>(p) = 3 + xp.get(0)*xp.get(0);
128 } else if (down.isInside(xp) == true) {
129 dw_p.add();
130 dw_p.last().get<0>() = p.getKey();
131 domain.getProp<1>(p) = 1 + xp.get(0)*xp.get(0);
132 } else if (left.isInside(xp) == true) {
133 l_p.add();
134 l_p.last().get<0>() = p.getKey();
135 domain.getProp<1>(p) = 1 + 2*xp.get(1)*xp.get(1);
136 } else if (right.isInside(xp) == true) {
137 r_p.add();
138 r_p.last().get<0>() = p.getKey();
139 domain.getProp<1>(p) = 2 + 2*xp.get(1)*xp.get(1);
140 } else {
141 bulk.add();
142 bulk.last().get<0>() = p.getKey();
143 }
144
145
146 ++it2;
147 }
148
149
150 auto eq1 = Lap(v);
151
152 Solver.impose(eq1, bulk, 6);
153 Solver.impose(v, up_p, RHS);
154 Solver.impose(v, dw_p, RHS);
155 Solver.impose(v, l_p, prop_id<1>());
156 Solver.impose(v, r_p, prop_id<1>());
157 Solver.solve(v);
158 anasol=Lap(v);
159
160 double worst1 = 0.0;
161
162 it2 = domain.getDomainIterator();
163
164 while (it2.isNext()) {
165 auto p = it2.get();
166 if (fabs(domain.getProp<0>(p) - domain.getProp<2>(p)) >= worst1) {
167 worst1 = fabs(domain.getProp<0>(p) - domain.getProp<2>(p));
168 }
169
170 domain.getProp<1>(p) = fabs(domain.getProp<0>(p) - domain.getProp<2>(p));
171
172 ++it2;
173 }
174
175 domain.write("particles");
176 BOOST_REQUIRE(worst1 < 0.03);
177 }
178
179
181 BOOST_AUTO_TEST_CASE(dcpse_poisson_Robin_anal) {
182// int rank;
183// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
184 const size_t sz[2] = {81,81};
185 Box<2, double> box({0, 0}, {0.5, 0.5});
186 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
187 double spacing = box.getHigh(0) / (sz[0] - 1);
188 Ghost<2, double> ghost(spacing * 3.1);
189 double rCut = 3.1 * spacing;
190 BOOST_TEST_MESSAGE("Init vector_dist...");
191
192 vector_dist_gpu<2, double, aggregate<double,double,double,double,double,double>> domain(0, box, bc, ghost);
193
194
195 //Init_DCPSE(domain)
196 BOOST_TEST_MESSAGE("Init domain...");
197
198 auto it = domain.getGridIterator(sz);
199 while (it.isNext()) {
200 domain.add();
201
202 auto key = it.get();
203 double x = key.get(0) * it.getSpacing(0);
204 domain.getLastPos()[0] = x;
205 double y = key.get(1) * it.getSpacing(1);
206 domain.getLastPos()[1] = y;
207
208 ++it;
209 }
210
211 // Add multi res patch 1
212
213 {
214 const size_t sz2[2] = {40,40};
215 Box<2,double> bx({0.25 + it.getSpacing(0)/4.0,0.25 + it.getSpacing(0)/4.0},{sz2[0]*it.getSpacing(0)/2.0 + 0.25 + it.getSpacing(0)/4.0, sz2[1]*it.getSpacing(0)/2.0 + 0.25 + it.getSpacing(0)/4.0});
216 openfpm::vector<size_t> rem;
217
218 auto it = domain.getDomainIterator();
219
220 while (it.isNext())
221 {
222 auto k = it.get();
223
224 Point<2,double> xp = domain.getPos(k);
225
226 if (bx.isInside(xp) == true)
227 {
228 rem.add(k.getKey());
229 }
230
231 ++it;
232 }
233
234 domain.remove(rem);
235
236 auto it2 = domain.getGridIterator(sz2);
237 while (it2.isNext()) {
238 domain.add();
239
240 auto key = it2.get();
241 double x = key.get(0) * spacing/2.0 + 0.25 + spacing/4.0;
242 domain.getLastPos()[0] = x;
243 double y = key.get(1) * spacing/2.0 + 0.25 + spacing/4.0;
244 domain.getLastPos()[1] = y;
245
246 ++it2;
247 }
248 }
249
250 // Add multi res patch 2
251
252 {
253 const size_t sz2[2] = {40,40};
254 Box<2,double> bx({0.25 + 21.0*spacing/8.0,0.25 + 21.0*spacing/8.0},{sz2[0]*spacing/4.0 + 0.25 + 21.0*spacing/8.0, sz2[1]*spacing/4.0 + 0.25 + 21*spacing/8.0});
255 openfpm::vector<size_t> rem;
256
257 auto it = domain.getDomainIterator();
258
259 while (it.isNext())
260 {
261 auto k = it.get();
262
263 Point<2,double> xp = domain.getPos(k);
264
265 if (bx.isInside(xp) == true)
266 {
267 rem.add(k.getKey());
268 }
269
270 ++it;
271 }
272
273 domain.remove(rem);
274
275 auto it2 = domain.getGridIterator(sz2);
276 while (it2.isNext()) {
277 domain.add();
278
279 auto key = it2.get();
280 double x = key.get(0) * spacing/4.0 + 0.25 + 21*spacing/8.0;
281 domain.getLastPos()[0] = x;
282 double y = key.get(1) * spacing/4.0 + 0.25 + 21*spacing/8.0;
283 domain.getLastPos()[1] = y;
284
285 ++it2;
286 }
287 }
288
290
291 BOOST_TEST_MESSAGE("Sync domain across processors...");
292
293 domain.map();
294 domain.ghost_get<0>();
295
296 Derivative_x_gpu Dx(domain, 2, rCut/3.0 ,1.9,support_options::N_PARTICLES);
297 Derivative_y_gpu Dy(domain, 2, rCut/3.0,1.9,support_options::N_PARTICLES);
298 Laplacian_gpu Lap(domain, 2, rCut/3.0 ,1.9,support_options::N_PARTICLES);
299
300 openfpm::vector<aggregate<int>> bulk;
301 openfpm::vector<aggregate<int>> up_p;
302 openfpm::vector<aggregate<int>> dw_p;
303 openfpm::vector<aggregate<int>> l_p;
304 openfpm::vector<aggregate<int>> r_p;
305 openfpm::vector<aggregate<int>> ref_p;
306
307 auto v = getV<0>(domain);
308 auto RHS=getV<1>(domain);
309 auto sol = getV<2>(domain);
310 auto anasol = getV<3>(domain);
311 auto err = getV<4>(domain);
312 auto DCPSE_sol=getV<5>(domain);
313
314 // Here fill me
315
316 Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
317 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
318
319 Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
320 {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
321
322 Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
323 {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
324
325 Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
326 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
327
328 openfpm::vector<Box<2, double>> boxes;
329 boxes.add(up);
330 boxes.add(down);
331 boxes.add(left);
332 boxes.add(right);
333
334 // Create a writer and write
335 VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
336 vtk_box.add(boxes);
337 //vtk_box.write("vtk_box.vtk");
338
339
340 auto it2 = domain.getDomainIterator();
341
342 while (it2.isNext()) {
343 auto p = it2.get();
344 Point<2, double> xp = domain.getPos(p);
345 if (up.isInside(xp) == true) {
346 up_p.add();
347 up_p.last().get<0>() = p.getKey();
348 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
349 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
350 } else if (down.isInside(xp) == true) {
351 dw_p.add();
352 dw_p.last().get<0>() = p.getKey();
353 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
354 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
355
356 } else if (left.isInside(xp) == true) {
357 l_p.add();
358 l_p.last().get<0>() = p.getKey();
359 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
360 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
361
362 } else if (right.isInside(xp) == true) {
363 r_p.add();
364 r_p.last().get<0>() = p.getKey();
365 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
366 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
367
368 } else {
369 bulk.add();
370 bulk.last().get<0>() = p.getKey();
371 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
372 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
373 }
374 ++it2;
375 }
376
377 domain.ghost_get<1,3>();
378
379 DCPSE_scheme_gpu<equations2d1_gpu,decltype(domain)> Solver( domain);
380 auto Poisson = Lap(v);
381 auto D_x = Dx(v);
382 auto D_y = Dy(v);
383 Solver.impose(Poisson, bulk, prop_id<1>());
384 Solver.impose(D_y, up_p, 0);
385 Solver.impose(D_x, r_p, 0);
386 Solver.impose(v, dw_p, 0);
387 Solver.impose(v, l_p, 0);
388
389 petsc_solver<double> solver;
390
391 solver.setPreconditioner(PCBJACOBI);
392 solver.setRestart(500);
393
394 Solver.solve_with_solver(solver,sol);
395
396 //solver.print_preconditioner();
397
398 domain.ghost_get<2>();
399
400 DCPSE_sol=Lap(sol);
401 domain.ghost_get<5>();
402
403 double worst1 = 0.0;
404
405 v=abs(DCPSE_sol-RHS);
406
407 for(int j=0;j<bulk.size();j++)
408 { auto p=bulk.get<0>(j);
409 if (fabs(domain.getProp<3>(p) - domain.getProp<2>(p)) >= worst1) {
410 worst1 = fabs(domain.getProp<3>(p) - domain.getProp<2>(p));
411 }
412 domain.getProp<4>(p) = fabs(domain.getProp<3>(p) - domain.getProp<2>(p));
413
414 }
415 //std::cout << "Maximum Analytic Error: " << worst1 << std::endl;
416
417 //domain.ghost_get<4>();
418 //domain.write("Robin_anasol");
419 BOOST_REQUIRE(worst1 < 0.03);
420
421 }
423 //81 0.00131586
424 //161 0.000328664
425 //320 8.30297e-05
426 //520 3.12398e-05
427 //1024 8.08087e-06
429 BOOST_AUTO_TEST_CASE(dcpse_poisson_Dirichlet_anal) {
430// int rank;
431// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
432 const size_t sz[2] = {81,81};
433 Box<2, double> box({0, 0}, {1, 1});
434 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
435 double spacing = box.getHigh(0) / (sz[0] - 1);
436 Ghost<2, double> ghost(spacing * 3.1);
437 double rCut = 3.1 * spacing;
438 BOOST_TEST_MESSAGE("Init vector_dist...");
439
440 vector_dist_gpu<2, double, aggregate<double,double,double,double,double,double>> domain(0, box, bc, ghost);
441
442
443 //Init_DCPSE(domain)
444 BOOST_TEST_MESSAGE("Init domain...");
445
446 auto it = domain.getGridIterator(sz);
447 while (it.isNext()) {
448 domain.add();
449
450 auto key = it.get();
451 double x = key.get(0) * it.getSpacing(0);
452 domain.getLastPos()[0] = x;
453 double y = key.get(1) * it.getSpacing(1);
454 domain.getLastPos()[1] = y;
455
456 ++it;
457 }
458 BOOST_TEST_MESSAGE("Sync domain across processors...");
459
460 domain.map();
461 domain.ghost_get<0>();
462
463 Derivative_x_gpu Dx(domain, 2, rCut,1.9,support_options::RADIUS);
464 Derivative_y_gpu Dy(domain, 2, rCut,1.9,support_options::RADIUS);
465 Laplacian_gpu Lap(domain, 2, rCut, 1.9,support_options::RADIUS);
466
467 openfpm::vector<aggregate<int>> bulk;
468 openfpm::vector<aggregate<int>> bulkF;
469 openfpm::vector<aggregate<int>> up_p;
470 openfpm::vector<aggregate<int>> dw_p;
471 openfpm::vector<aggregate<int>> l_p;
472 openfpm::vector<aggregate<int>> r_p;
473 openfpm::vector<aggregate<int>> ref_p;
474
475 auto v = getV<0>(domain);
476 auto RHS=getV<1>(domain);
477 auto sol = getV<2>(domain);
478 auto anasol = getV<3>(domain);
479 auto err = getV<4>(domain);
480 auto DCPSE_sol=getV<5>(domain);
481
482 // Here fill me
483
484 Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
485 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
486
487 Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
488 {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
489
490 Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
491 {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
492
493 Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
494 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
495
496 openfpm::vector<Box<2, double>> boxes;
497 boxes.add(up);
498 boxes.add(down);
499 boxes.add(left);
500 boxes.add(right);
501
502 // Create a writer and write
503 VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
504 vtk_box.add(boxes);
505 //vtk_box.write("vtk_box.vtk");
506
507
508 auto it2 = domain.getDomainIterator();
509
510 while (it2.isNext()) {
511 auto p = it2.get();
512 Point<2, double> xp = domain.getPos(p);
513 if (up.isInside(xp) == true) {
514 up_p.add();
515 up_p.last().get<0>() = p.getKey();
516 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
517 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
518 bulkF.add();
519 bulkF.last().get<0>() = p.getKey();
520 } else if (down.isInside(xp) == true) {
521 dw_p.add();
522 dw_p.last().get<0>() = p.getKey();
523 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
524 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
525 bulkF.add();
526 bulkF.last().get<0>() = p.getKey();
527
528 } else if (left.isInside(xp) == true) {
529 l_p.add();
530 l_p.last().get<0>() = p.getKey();
531 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
532 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
533 bulkF.add();
534 bulkF.last().get<0>() = p.getKey();
535
536 } else if (right.isInside(xp) == true) {
537 r_p.add();
538 r_p.last().get<0>() = p.getKey();
539 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
540 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
541 bulkF.add();
542 bulkF.last().get<0>() = p.getKey();
543
544 } else {
545 bulk.add();
546 bulk.last().get<0>() = p.getKey();
547 domain.getProp<1>(p) = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
548 domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
549 bulkF.add();
550 bulkF.last().get<0>() = p.getKey();
551 }
552 ++it2;
553 }
554 DCPSE_scheme_gpu<equations2d1_gpu,decltype(domain)> Solver( domain);
555 auto Poisson = Lap(v);
556 auto D_x = Dx(v);
557 auto D_y = Dy(v);
558 Solver.impose(Poisson, bulk, prop_id<1>());
559 Solver.impose(v, up_p, prop_id<1>());
560 Solver.impose(v, r_p, prop_id<1>());
561 Solver.impose(v, dw_p, prop_id<1>());
562 Solver.impose(v, l_p, prop_id<1>());
563 Solver.solve(sol);
564 DCPSE_sol=Lap(sol);
565
566
567 for (int j = 0; j < up_p.size(); j++) {
568 auto p = up_p.get<0>(j);
569 domain.getProp<5>(p) = 0;
570 }
571 for (int j = 0; j < dw_p.size(); j++) {
572 auto p = dw_p.get<0>(j);
573 domain.getProp<5>(p) = 0;
574 }
575 for (int j = 0; j < l_p.size(); j++) {
576 auto p = l_p.get<0>(j);
577 domain.getProp<5>(p) = 0;
578 }
579 for (int j = 0; j < r_p.size(); j++) {
580 auto p = r_p.get<0>(j);
581 domain.getProp<5>(p) = 0;
582 }
583
584 double worst1 = 0.0;
585
586 v=abs(DCPSE_sol-RHS);
587
588 for(int j=0;j<bulkF.size();j++)
589 { auto p=bulkF.get<0>(j);
590 if (fabs(domain.getProp<3>(p) - domain.getProp<2>(p)) >= worst1) {
591 worst1 = fabs(domain.getProp<3>(p) - domain.getProp<2>(p));
592 }
593 domain.getProp<4>(p) = fabs(domain.getProp<3>(p) - domain.getProp<2>(p));
594
595 }
596 // std::cout << "Maximum Analytic Error: " << worst1 << std::endl;
597
598 BOOST_REQUIRE(worst1 < 0.03);
599
600 // domain.write("Dirichlet_anasol");
601 }
603
604 BOOST_AUTO_TEST_CASE(dcpse_poisson_Periodic) {
605 //https://fenicsproject.org/docs/dolfin/1.4.0/python/demo/documented/periodic/python/documentation.html
606 // int rank;
607 // MPI_Comm_rank(MPI_COMM_WORLD, &rank);
608 const size_t sz[2] = {31,31};
609 Box<2, double> box({0, 0}, {1.0, 1.0});
610 size_t bc[2] = {PERIODIC, NON_PERIODIC};
611 double spacing = box.getHigh(0) / (sz[0] - 1);
612 Ghost<2, double> ghost(spacing * 3.1);
613 double rCut = 3.1 * spacing;
614 BOOST_TEST_MESSAGE("Init vector_dist...");
615
616 vector_dist_gpu<2, double, aggregate<double,double,double,double,double,VectorS<2, double>>> domain(0, box, bc, ghost);
617
618 //Init_DCPSE(domain)
619 BOOST_TEST_MESSAGE("Init domain...");
620
621 auto it = domain.getGridIterator(sz);
622 while (it.isNext()) {
623 domain.add();
624
625 auto key = it.get();
626 double x = key.get(0) * it.getSpacing(0);
627 domain.getLastPos()[0] = x;
628 double y = key.get(1) * it.getSpacing(1)*0.99999;
629 domain.getLastPos()[1] = y;
630
631 ++it;
632 }
633 BOOST_TEST_MESSAGE("Sync domain across processors...");
634
635 domain.map();
636 domain.ghost_get<0>();
637
638
639 Laplacian_gpu Lap(domain, 2, rCut, 1.9, support_options::RADIUS);
640
641 DCPSE_scheme_gpu<equations2d1p_gpu,decltype(domain)> Solver( domain);
642
643 openfpm::vector<aggregate<int>> bulk;
644 openfpm::vector<aggregate<int>> up_p;
645 openfpm::vector<aggregate<int>> dw_p;
646 openfpm::vector<aggregate<int>> l_p;
647 openfpm::vector<aggregate<int>> r_p;
648
649 auto v = getV<0>(domain);
650 auto RHS=getV<1>(domain);
651 auto sol = getV<2>(domain);
652 auto anasol = getV<3>(domain);
653 auto err = getV<4>(domain);
654 auto u = getV<5>(domain);
655
656 // Here fill me
657
658 Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
659 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
660
661 Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
662 {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
663
664 Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
665 {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
666
667 Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
668 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
669
670 openfpm::vector<Box<2, double>> boxes;
671 boxes.add(up);
672 boxes.add(down);
673 boxes.add(left);
674 boxes.add(right);
675
676 // Create a writer and write
677 VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
678 vtk_box.add(boxes);
679 //vtk_box.write("vtk_box.vtk");
680
681
682 auto it2 = domain.getDomainIterator();
683
684 while (it2.isNext()) {
685 auto p = it2.get();
686 Point<2, double> xp = domain.getPos(p);
687 //domain.getProp<3>(p)=1+xp[0]*xp[0]+2*xp[1]*xp[1];
688 if (up.isInside(xp) == true) {
689 up_p.add();
690 up_p.last().get<0>() = p.getKey();
691 domain.getProp<1>(p) = 0;
692 } else if (down.isInside(xp) == true) {
693 dw_p.add();
694 dw_p.last().get<0>() = p.getKey();
695 domain.getProp<1>(p) = 0;
696 } else {
697 bulk.add();
698 bulk.last().get<0>() = p.getKey();
699 domain.getProp<1>(p) = xp.get(0)*sin(5*M_PI*xp.get(1))+exp(-((xp.get(0)-0.5)*(xp.get(0)-0.5)+(xp.get(1)-0.5)*(xp.get(1)-0.5))/0.02);
700 }
701
702 ++it2;
703 }
704
705 domain.ghost_get<1>();
706 auto Poisson = -Lap(v);
707 Solver.impose(Poisson, bulk, prop_id<1>());
708 Solver.impose(v, up_p, 0);
709 Solver.impose(v, dw_p, 0);
710 Solver.solve(v);
711
712 domain.ghost_get<0>();
713 anasol=-Lap(v);
714 double worst1 = 0.0;
715
716 for(int j=0;j<bulk.size();j++)
717 { auto p=bulk.get<0>(j);
718 if (fabs(domain.getProp<3>(p) - domain.getProp<1>(p)) >= worst1) {
719 worst1 = fabs(domain.getProp<3>(p) - domain.getProp<1>(p));
720 }
721 domain.getProp<4>(p) = fabs(domain.getProp<1>(p) - domain.getProp<3>(p));
722
723 }
724 //Auto Error
725 BOOST_REQUIRE(worst1 < 1.0);
726
727 //domain.write("Poisson_Periodic");
728 }
730
731 BOOST_AUTO_TEST_CASE(dcpse_poisson_Robin) {
732 //http://e6.ijs.si/medusa/wiki/index.php/Poisson%27s_equation#Full_Neumann_boundary_conditions
733// int rank;
734// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
735 const size_t sz[2] = {31,31};
736 Box<2, double> box({0, 0}, {1.0, 1.0});
737 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
738 double spacing = box.getHigh(0) / (sz[0] - 1);
739 Ghost<2, double> ghost(spacing * 3);
740 double rCut = 2.0 * spacing;
741 BOOST_TEST_MESSAGE("Init vector_dist...");
742
743 vector_dist_gpu<2, double, aggregate<double,double,double,double,double,VectorS<2, double>>> domain(0, box, bc, ghost);
744
745
746 //Init_DCPSE(domain)
747 BOOST_TEST_MESSAGE("Init domain...");
748
749 auto it = domain.getGridIterator(sz);
750 while (it.isNext()) {
751 domain.add();
752
753 auto key = it.get();
754 double x = key.get(0) * it.getSpacing(0);
755 domain.getLastPos()[0] = x;
756 double y = key.get(1) * it.getSpacing(1);
757 domain.getLastPos()[1] = y;
758
759 ++it;
760 }
761 BOOST_TEST_MESSAGE("Sync domain across processors...");
762
763 domain.map();
764 domain.ghost_get<0>();
765
766 Derivative_y_gpu Dy(domain, 2, rCut,2,support_options::N_PARTICLES);
767 Laplacian_gpu Lap(domain, 2, rCut, 3,support_options::N_PARTICLES);
768
769 DCPSE_scheme_gpu<equations2d1_gpu,decltype(domain)> Solver(domain);
770
771 openfpm::vector<aggregate<int>> bulk;
772 openfpm::vector<aggregate<int>> up_p;
773 openfpm::vector<aggregate<int>> dw_p;
774 openfpm::vector<aggregate<int>> l_p;
775 openfpm::vector<aggregate<int>> r_p;
776
777 auto v = getV<0>(domain);
778 auto sol = getV<2>(domain);
779 auto anasol = getV<3>(domain);
780 auto err = getV<4>(domain);
781 auto u = getV<5>(domain);
782
783 // Here fill me
784
785 Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
786 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
787
788 Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
789 {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
790
791 Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
792 {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
793
794 Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
795 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
796
797 openfpm::vector<Box<2, double>> boxes;
798 boxes.add(up);
799 boxes.add(down);
800 boxes.add(left);
801 boxes.add(right);
802
803 // Create a writer and write
804 VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
805 vtk_box.add(boxes);
806 //vtk_box.write("vtk_box.vtk");
807
808
809 auto it2 = domain.getDomainIterator();
810
811 while (it2.isNext()) {
812 auto p = it2.get();
813 Point<2, double> xp = domain.getPos(p);
814 //domain.getProp<3>(p)=1+xp[0]*xp[0]+2*xp[1]*xp[1];
815 if (up.isInside(xp) == true) {
816 up_p.add();
817 up_p.last().get<0>() = p.getKey();
818 domain.getProp<1>(p) = sin(5.0*xp.get(0));
819 } else if (down.isInside(xp) == true) {
820 dw_p.add();
821 dw_p.last().get<0>() = p.getKey();
822 domain.getProp<1>(p) = sin(5.0*xp.get(0));
823 } else if (left.isInside(xp) == true) {
824 l_p.add();
825 l_p.last().get<0>() = p.getKey();
826 domain.getProp<1>(p) = sin(5.0*xp.get(0));
827 } else if (right.isInside(xp) == true) {
828 r_p.add();
829 r_p.last().get<0>() = p.getKey();
830 domain.getProp<1>(p) = sin(5.0*xp.get(0));
831 } else {
832 bulk.add();
833 bulk.last().get<0>() = p.getKey();
834 domain.getProp<1>(p) = -10.0*exp(-((xp.get(0)-0.5)*(xp.get(0)-0.5)+(xp.get(1)-0.5)*(xp.get(1)-0.5))/0.02);
835 }
836
837 ++it2;
838 }
839
840 petsc_solver<double> pet_sol;
841 pet_sol.setPreconditioner(PCNONE);
842
843 auto Poisson = Lap(v);
844 auto D_y = Dy(v);
845 Solver.impose(Poisson, bulk, prop_id<1>());
846 Solver.impose(D_y, up_p, prop_id<1>());
847 Solver.impose(-D_y, dw_p, prop_id<1>());
848 Solver.impose(v, l_p, 0);
849 Solver.impose(v, r_p, 0);
850
851 Solver.solve_with_solver(pet_sol,sol);
852 domain.ghost_get<2>();
853
854 anasol=Lap(sol);
855 double worst1 = 0.0;
856
857 for(int j=0;j<bulk.size();j++)
858 { auto p=bulk.get<0>(j);
859 if (fabs(domain.getProp<3>(p) - domain.getProp<1>(p)) >= worst1) {
860 worst1 = fabs(domain.getProp<3>(p) - domain.getProp<1>(p));
861 }
862 domain.getProp<4>(p) = fabs(domain.getProp<1>(p) - domain.getProp<3>(p));
863
864 }
865 //Auto Error
866 BOOST_REQUIRE(worst1 < 1.0);
867
868 //std::cout << "WORST: " << worst1 << std::endl;
869
870 //domain.write("Mixed");
871 }
873
875
876
877 BOOST_AUTO_TEST_CASE(dcpse_poisson_Neumann) {
878 //https://fenicsproject.org/docs/dolfin/1.4.0/python/demo/documented/neumann-poisson/python/documentation.html
879// int rank;
880// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
881 const size_t sz[2] = {31,31};
882 Box<2, double> box({0, 0}, {1.0, 1.0});
883 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
884 double spacing = box.getHigh(0) / (sz[0] - 1);
885 double rCut = 3.1 * spacing;
886 Ghost<2, double> ghost(spacing * 3.1);
887 BOOST_TEST_MESSAGE("Init vector_dist...");
888
889 vector_dist_gpu<2, double, aggregate<double,double,double,double,double>> domain(0, box, bc, ghost);
890
891
892 //Init_DCPSE(domain)
893 BOOST_TEST_MESSAGE("Init domain...");
894
895 auto it = domain.getGridIterator(sz);
896 while (it.isNext()) {
897
898 domain.add();
899
900 auto key = it.get();
901 double x = key.get(0) * it.getSpacing(0);
902 domain.getLastPos()[0] = x;
903 double y = key.get(1) * it.getSpacing(1);
904 domain.getLastPos()[1] = y;
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_gpu Dx(domain, 2, rCut,1.9,support_options::N_PARTICLES);
914 Derivative_y_gpu Dy(domain, 2, rCut,1.9,support_options::N_PARTICLES);
915 Laplacian_gpu Lap(domain, 2, rCut, 1.9,support_options::N_PARTICLES);
916
917 petsc_solver<double> solver;
918 solver.setRestart(500);
919 solver.setSolver(KSPGMRES);
920 solver.setPreconditioner(PCSVD);
921
922 openfpm::vector<aggregate<int>> bulk;
923 openfpm::vector<aggregate<int>> up_p;
924 openfpm::vector<aggregate<int>> dw_p;
925 openfpm::vector<aggregate<int>> l_p;
926 openfpm::vector<aggregate<int>> r_p;
927
928 auto v = getV<0>(domain);
929 //auto RHS=getV<1>(domain);
930 auto sol = getV<2>(domain);
931 auto anasol = getV<3>(domain);
932 auto err = getV<4>(domain);
933
934 // Here fill me
935
936 Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
937 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
938
939 Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
940 {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
941
942 Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
943 {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
944
945 Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
946 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
947
948 openfpm::vector<Box<2, double>> boxes;
949 boxes.add(up);
950 boxes.add(down);
951 boxes.add(left);
952 boxes.add(right);
953
954 // Create a writer and write
955 VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
956 vtk_box.add(boxes);
957 //vtk_box.write("vtk_box.vtk");
958
959
960 auto it2 = domain.getDomainIterator();
961
962 while (it2.isNext()) {
963 auto p = it2.get();
964 Point<2, double> xp = domain.getPos(p);
965 //domain.getProp<3>(p)=1+xp[0]*xp[0]+2*xp[1]*xp[1];
966 if (up.isInside(xp) == true) {
967 up_p.add();
968 up_p.last().get<0>() = p.getKey();
969 domain.getProp<1>(p) = sin(5*xp.get(0));
970 } else if (down.isInside(xp) == true) {
971 dw_p.add();
972 dw_p.last().get<0>() = p.getKey();
973 domain.getProp<1>(p) = sin(5*xp.get(0));
974 } else if (left.isInside(xp) == true) {
975 l_p.add();
976 l_p.last().get<0>() = p.getKey();
977 domain.getProp<1>(p) = sin(5*xp.get(0));
978 } else if (right.isInside(xp) == true) {
979 r_p.add();
980 r_p.last().get<0>() = p.getKey();
981 domain.getProp<1>(p) = sin(5*xp.get(0));
982 } else {
983 bulk.add();
984 bulk.last().get<0>() = p.getKey();
985 domain.getProp<1>(p) = -10*exp(-((xp.get(0)-0.5)*(xp.get(0)-0.5)+(xp.get(1)-0.5)*(xp.get(1)-0.5))/0.02);
986 }
987
988 ++it2;
989 }
990
991 DCPSE_scheme_gpu<equations2d1_gpu,decltype(domain)> Solver(domain,options_solver::LAGRANGE_MULTIPLIER);
992 auto Poisson = -Lap(v);
993 auto D_x = Dx(v);
994 auto D_y = Dy(v);
995 Solver.impose(Poisson, bulk, prop_id<1>());
996 Solver.impose(D_y, up_p, prop_id<1>());
997 Solver.impose(-D_y, dw_p, prop_id<1>());
998 Solver.impose(-D_x, l_p, prop_id<1>());
999 Solver.impose(D_x, r_p, prop_id<1>());
1000 Solver.solve_with_solver(solver,sol);
1001
1002// Solver.solve(sol);
1003 domain.ghost_get<2>();
1004 anasol=-Lap(sol);
1005 double worst1 = 0.0;
1006
1007 for(int j=0;j<bulk.size();j++)
1008 { auto p=bulk.get<0>(j);
1009 if (fabs(domain.getProp<3>(p) - domain.getProp<1>(p)) >= worst1) {
1010 worst1 = fabs(domain.getProp<3>(p) - domain.getProp<1>(p));
1011 }
1012 domain.getProp<4>(p) = fabs(domain.getProp<1>(p) - domain.getProp<3>(p));
1013
1014 }
1015 //Auto Error
1016 BOOST_REQUIRE(worst1 < 1.0);
1017
1018 //domain.write("Neumann");
1019 }
1020
1021
1022 BOOST_AUTO_TEST_CASE(dcpse_slice_solver) {
1023// int rank;
1024// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1025 const size_t sz[2] = {31,31};
1026 Box<2, double> box({0, 0}, {1, 1});
1027 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
1028 double spacing = box.getHigh(0) / (sz[0] - 1);
1029 double rCut = 3.1 * spacing;
1030 Ghost<2, double> ghost(rCut);
1031 BOOST_TEST_MESSAGE("Init vector_dist...");
1032
1033 vector_dist_gpu<2, double, aggregate<VectorS<2, double>,VectorS<2, double>,VectorS<2, double>,VectorS<2, double>,VectorS<2, double>,VectorS<2, double>>> domain(0, box, bc, ghost);
1034
1035
1036 //Init_DCPSE(domain)
1037 BOOST_TEST_MESSAGE("Init domain...");
1038
1039 auto it = domain.getGridIterator(sz);
1040 while (it.isNext()) {
1041 domain.add();
1042
1043 auto key = it.get();
1044 double x = key.get(0) * it.getSpacing(0);
1045 domain.getLastPos()[0] = x;
1046 double y = key.get(1) * it.getSpacing(1);
1047 domain.getLastPos()[1] = y;
1048
1049 ++it;
1050 }
1051 BOOST_TEST_MESSAGE("Sync domain across processors...");
1052
1053 domain.map();
1054 domain.ghost_get<0>();
1055
1056 Derivative_x Dx(domain, 2, rCut,1.9,support_options::RADIUS);
1057 Derivative_y Dy(domain, 2, rCut,1.9,support_options::RADIUS);
1058 Laplacian Lap(domain, 2, rCut,1.9,support_options::RADIUS);
1059
1060 openfpm::vector<aggregate<int>> bulk;
1061 openfpm::vector<aggregate<int>> up_p;
1062 openfpm::vector<aggregate<int>> dw_p;
1063 openfpm::vector<aggregate<int>> l_p;
1064 openfpm::vector<aggregate<int>> r_p;
1065 openfpm::vector<aggregate<int>> ref_p;
1066
1067 auto v = getV<0>(domain);
1068 auto RHS=getV<1>(domain);
1069 auto sol = getV<2>(domain);
1070 auto anasol = getV<3>(domain);
1071 auto err = getV<4>(domain);
1072 auto DCPSE_sol=getV<5>(domain);
1073
1074 // Here fill me
1075
1076 Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
1077 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
1078
1079 Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
1080 {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
1081
1082 Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
1083 {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
1084
1085 Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
1086 {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
1087
1088 openfpm::vector<Box<2, double>> boxes;
1089 boxes.add(up);
1090 boxes.add(down);
1091 boxes.add(left);
1092 boxes.add(right);
1093
1094 // Create a writer and write
1095 VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
1096 vtk_box.add(boxes);
1097 //vtk_box.write("vtk_box.vtk");
1098 // domain.write("Slice_anasol");
1099
1100
1101 auto it2 = domain.getDomainIterator();
1102
1103 while (it2.isNext()) {
1104 auto p = it2.get();
1105 Point<2, double> xp = domain.getPos(p);
1106 if (up.isInside(xp) == true) {
1107 up_p.add();
1108 up_p.last().get<0>() = p.getKey();
1109 domain.getProp<1>(p)[0] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1110 domain.getProp<3>(p)[0] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1111
1112 domain.getProp<1>(p)[1] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1113 domain.getProp<3>(p)[1] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1114
1115
1116 } else if (down.isInside(xp) == true) {
1117 dw_p.add();
1118 dw_p.last().get<0>() = p.getKey();
1119 domain.getProp<1>(p)[0] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1120 domain.getProp<3>(p)[0] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1121
1122 domain.getProp<1>(p)[1] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1123 domain.getProp<3>(p)[1] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1124
1125 } else if (left.isInside(xp) == true) {
1126 l_p.add();
1127 l_p.last().get<0>() = p.getKey();
1128 domain.getProp<1>(p)[0] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1129 domain.getProp<3>(p)[0] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1130
1131 domain.getProp<1>(p)[1] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1132 domain.getProp<3>(p)[1] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1133
1134 } else if (right.isInside(xp) == true) {
1135 r_p.add();
1136 r_p.last().get<0>() = p.getKey();
1137 domain.getProp<1>(p)[0] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1138 domain.getProp<3>(p)[0] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1139
1140 domain.getProp<1>(p)[1] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1141 domain.getProp<3>(p)[1] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1142
1143 } else {
1144 bulk.add();
1145 bulk.last().get<0>() = p.getKey();
1146 domain.getProp<1>(p)[0] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1147 domain.getProp<3>(p)[0] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1148
1149 domain.getProp<1>(p)[1] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1150 domain.getProp<3>(p)[1] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1151 }
1152 ++it2;
1153 }
1154
1155 eq_id vx,vy;
1156
1157 vx.setId(0);
1158 vy.setId(1);
1159
1160 DCPSE_scheme_gpu<equations2d2_gpu,decltype(domain)> Solver( domain);
1161 auto Poisson0 = Lap(v[0]);
1162 auto Poisson1 = Lap(v[1]);
1163 //auto D_x = Dx(v[1]);
1164 //auto D_y = Dy(v[1]);
1165 Solver.impose(Poisson0, bulk, RHS[0],vx);
1166 Solver.impose(Poisson1, bulk, RHS[1],vy);
1167 Solver.impose(v[0], up_p, RHS[0],vx);
1168 Solver.impose(v[1], up_p, RHS[1],vy);
1169 Solver.impose(v[0], r_p, RHS[0],vx);
1170 Solver.impose(v[1], r_p, RHS[1],vy);
1171 Solver.impose(v[0], dw_p, RHS[0],vx);
1172 Solver.impose(v[1], dw_p, RHS[1],vy);
1173 Solver.impose(v[0], l_p, RHS[0],vx);
1174 Solver.impose(v[1], l_p, RHS[1],vy);
1175 Solver.solve(sol[0],sol[1]);
1176 DCPSE_sol=Lap(sol);
1177 double worst1 = 0.0;
1178 double worst2 = 0.0;
1179
1180
1181 v=sol-RHS;
1182
1183 for(int j=0;j<bulk.size();j++)
1184 { auto p=bulk.get<0>(j);
1185 if (fabs(domain.getProp<3>(p)[0] - domain.getProp<2>(p)[0]) >= worst1) {
1186 worst1 = fabs(domain.getProp<3>(p)[0] - domain.getProp<2>(p)[0]);
1187 }
1188 domain.getProp<4>(p)[0] = fabs(domain.getProp<3>(p)[0] - domain.getProp<2>(p)[0]);
1189
1190 }
1191 for(int j=0;j<bulk.size();j++)
1192 { auto p=bulk.get<0>(j);
1193 if (fabs(domain.getProp<3>(p)[1] - domain.getProp<2>(p)[1]) >= worst2) {
1194 worst2 = fabs(domain.getProp<3>(p)[1] - domain.getProp<2>(p)[1]);
1195 }
1196 domain.getProp<4>(p)[1] = fabs(domain.getProp<3>(p)[1] - domain.getProp<2>(p)[1]);
1197
1198 }
1199 //std::cout << "Maximum Analytic Error in slice x: " << worst1 << std::endl;
1200 //std::cout << "Maximum Analytic Error in slice y: " << worst2 << std::endl;
1201 //domain.write("Slice_anasol");
1202 BOOST_REQUIRE(worst1 < 0.03);
1203 BOOST_REQUIRE(worst2 < 0.03);
1204
1205 }
1206
1207
1208BOOST_AUTO_TEST_SUITE_END()
1209
1210
Box
This class represent an N-dimensional box.
Definition Box.hpp:61
Lap
Laplacian second order on h (spacing)
Definition Laplacian.hpp:23
Point
This class implement the point shape in an N-dimensional space.
Definition Point.hpp:28
Point::get
__device__ __host__ const T & get(unsigned int i) const
Get coordinate.
Definition Point.hpp:172
eq_id
Definition eq_solve_common.hpp:15
openfpm::vector
Implementation of 1-D std::vector like structure.
Definition map_vector.hpp:203
openfpm::vector::size
size_t size()
Stub size.
Definition map_vector.hpp:211
petsc_solver
In case T does not match the PETSC precision compilation create a stub structure.
Definition petsc_solver.hpp:56
vector_dist
Distributed vector.
Definition vector_dist.hpp:305
prop_id
Definition eq_solve_common.hpp:12