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