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