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