OpenFPM_pdata  4.1.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_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 
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 * 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  Solver.solve_with_solver(solver,sol);
1005 
1006 // Solver.solve(sol);
1007  domain.ghost_get<2>();
1008  anasol=-Lap(sol);
1009  double worst1 = 0.0;
1010 
1011  for(int j=0;j<bulk.size();j++)
1012  { auto p=bulk.get<0>(j);
1013  if (fabs(domain.getProp<3>(p) - domain.getProp<1>(p)) >= worst1) {
1014  worst1 = fabs(domain.getProp<3>(p) - domain.getProp<1>(p));
1015  }
1016  domain.getProp<4>(p) = fabs(domain.getProp<1>(p) - domain.getProp<3>(p));
1017 
1018  }
1019  //Auto Error
1020  BOOST_REQUIRE(worst1 < 1.0);
1021 
1022  //domain.write("Neumann");
1023  }
1024 
1025 
1026  BOOST_AUTO_TEST_CASE(dcpse_slice_solver) {
1027 // int rank;
1028 // MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1029  const size_t sz[2] = {31,31};
1030  Box<2, double> box({0, 0}, {1, 1});
1031  size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
1032  double spacing = box.getHigh(0) / (sz[0] - 1);
1033  double rCut = 3.1 * spacing;
1034  Ghost<2, double> ghost(rCut);
1035  BOOST_TEST_MESSAGE("Init vector_dist...");
1036 
1037  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);
1038 
1039 
1040  //Init_DCPSE(domain)
1041  BOOST_TEST_MESSAGE("Init domain...");
1042 
1043  auto it = domain.getGridIterator(sz);
1044  while (it.isNext()) {
1045  domain.add();
1046 
1047  auto key = it.get();
1048  double x = key.get(0) * it.getSpacing(0);
1049  domain.getLastPos()[0] = x;
1050  double y = key.get(1) * it.getSpacing(1);
1051  domain.getLastPos()[1] = y;
1052 
1053  ++it;
1054  }
1055  BOOST_TEST_MESSAGE("Sync domain across processors...");
1056 
1057  domain.map();
1058  domain.ghost_get<0>();
1059 
1060  Derivative_x Dx(domain, 2, rCut,1.9,support_options::RADIUS);
1061  Derivative_y Dy(domain, 2, rCut,1.9,support_options::RADIUS);
1062  Laplacian Lap(domain, 2, rCut,1.9,support_options::RADIUS);
1063 
1064 
1065  openfpm::vector<aggregate<int>> bulk;
1066  openfpm::vector<aggregate<int>> up_p;
1067  openfpm::vector<aggregate<int>> dw_p;
1068  openfpm::vector<aggregate<int>> l_p;
1069  openfpm::vector<aggregate<int>> r_p;
1070  openfpm::vector<aggregate<int>> ref_p;
1071 
1072  auto v = getV<0>(domain);
1073  auto RHS=getV<1>(domain);
1074  auto sol = getV<2>(domain);
1075  auto anasol = getV<3>(domain);
1076  auto err = getV<4>(domain);
1077  auto DCPSE_sol=getV<5>(domain);
1078 
1079  // Here fill me
1080 
1081  Box<2, double> up({box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0},
1082  {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0});
1083 
1084  Box<2, double> down({box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0},
1085  {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0});
1086 
1087  Box<2, double> left({box.getLow(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
1088  {box.getLow(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
1089 
1090  Box<2, double> right({box.getHigh(0) - spacing / 2.0, box.getLow(1) + spacing / 2.0},
1091  {box.getHigh(0) + spacing / 2.0, box.getHigh(1) - spacing / 2.0});
1092 
1093  openfpm::vector<Box<2, double>> boxes;
1094  boxes.add(up);
1095  boxes.add(down);
1096  boxes.add(left);
1097  boxes.add(right);
1098 
1099  // Create a writer and write
1100  VTKWriter<openfpm::vector<Box<2, double>>, VECTOR_BOX> vtk_box;
1101  vtk_box.add(boxes);
1102  //vtk_box.write("vtk_box.vtk");
1103  // domain.write("Slice_anasol");
1104 
1105 
1106  auto it2 = domain.getDomainIterator();
1107 
1108  while (it2.isNext()) {
1109  auto p = it2.get();
1110  Point<2, double> xp = domain.getPos(p);
1111  if (up.isInside(xp) == true) {
1112  up_p.add();
1113  up_p.last().get<0>() = p.getKey();
1114  domain.getProp<1>(p)[0] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1115  domain.getProp<3>(p)[0] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1116 
1117  domain.getProp<1>(p)[1] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1118  domain.getProp<3>(p)[1] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1119 
1120 
1121  } else if (down.isInside(xp) == true) {
1122  dw_p.add();
1123  dw_p.last().get<0>() = p.getKey();
1124  domain.getProp<1>(p)[0] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1125  domain.getProp<3>(p)[0] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1126 
1127  domain.getProp<1>(p)[1] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1128  domain.getProp<3>(p)[1] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1129 
1130  } else if (left.isInside(xp) == true) {
1131  l_p.add();
1132  l_p.last().get<0>() = p.getKey();
1133  domain.getProp<1>(p)[0] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1134  domain.getProp<3>(p)[0] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1135 
1136  domain.getProp<1>(p)[1] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1137  domain.getProp<3>(p)[1] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1138 
1139  } else if (right.isInside(xp) == true) {
1140  r_p.add();
1141  r_p.last().get<0>() = p.getKey();
1142  domain.getProp<1>(p)[0] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1143  domain.getProp<3>(p)[0] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1144 
1145  domain.getProp<1>(p)[1] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1146  domain.getProp<3>(p)[1] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1147 
1148  } else {
1149  bulk.add();
1150  bulk.last().get<0>() = p.getKey();
1151  domain.getProp<1>(p)[0] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1152  domain.getProp<3>(p)[0] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1153 
1154  domain.getProp<1>(p)[1] = -2*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1155  domain.getProp<3>(p)[1] = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1));
1156  }
1157  ++it2;
1158  }
1159 
1160  eq_id vx,vy;
1161 
1162  vx.setId(0);
1163  vy.setId(1);
1164 
1165  DCPSE_scheme<equations2d2,decltype(domain)> Solver( domain);
1166  auto Poisson0 = Lap(v[0]);
1167  auto Poisson1 = Lap(v[1]);
1168  //auto D_x = Dx(v[1]);
1169  //auto D_y = Dy(v[1]);
1170  Solver.impose(Poisson0, bulk, RHS[0],vx);
1171  Solver.impose(Poisson1, bulk, RHS[1],vy);
1172  Solver.impose(v[0], up_p, RHS[0],vx);
1173  Solver.impose(v[1], up_p, RHS[1],vy);
1174  Solver.impose(v[0], r_p, RHS[0],vx);
1175  Solver.impose(v[1], r_p, RHS[1],vy);
1176  Solver.impose(v[0], dw_p, RHS[0],vx);
1177  Solver.impose(v[1], dw_p, RHS[1],vy);
1178  Solver.impose(v[0], l_p, RHS[0],vx);
1179  Solver.impose(v[1], l_p, RHS[1],vy);
1180  Solver.solve(sol[0],sol[1]);
1181  DCPSE_sol=Lap(sol);
1182  double worst1 = 0.0;
1183  double worst2 = 0.0;
1184 
1185 
1186  v=sol-RHS;
1187 
1188  for(int j=0;j<bulk.size();j++)
1189  { auto p=bulk.get<0>(j);
1190  if (fabs(domain.getProp<3>(p)[0] - domain.getProp<2>(p)[0]) >= worst1) {
1191  worst1 = fabs(domain.getProp<3>(p)[0] - domain.getProp<2>(p)[0]);
1192  }
1193  domain.getProp<4>(p)[0] = fabs(domain.getProp<3>(p)[0] - domain.getProp<2>(p)[0]);
1194 
1195  }
1196  for(int j=0;j<bulk.size();j++)
1197  { auto p=bulk.get<0>(j);
1198  if (fabs(domain.getProp<3>(p)[1] - domain.getProp<2>(p)[1]) >= worst2) {
1199  worst2 = fabs(domain.getProp<3>(p)[1] - domain.getProp<2>(p)[1]);
1200  }
1201  domain.getProp<4>(p)[1] = fabs(domain.getProp<3>(p)[1] - domain.getProp<2>(p)[1]);
1202 
1203  }
1204  //std::cout << "Maximum Analytic Error in slice x: " << worst1 << std::endl;
1205  //std::cout << "Maximum Analytic Error in slice y: " << worst2 << std::endl;
1206  //domain.write("Slice_anasol");
1207  BOOST_REQUIRE(worst1 < 0.03);
1208  BOOST_REQUIRE(worst2 < 0.03);
1209 
1210  }
1211 
1212 
1213 BOOST_AUTO_TEST_SUITE_END()
1214 #endif
1215 #endif
1216 
petsc_solver< double >::setPreconditioner
void setPreconditioner(PCType type)
Set the preconditioner of the linear solver.
Definition: petsc_solver.hpp:1046
equations2d1
Specify the general characteristic of system to solve.
Definition: EqnsStruct.hpp:13
Point
This class implement the point shape in an N-dimensional space.
Definition: Point.hpp:27
openfpm::vector::size
size_t size()
Stub size.
Definition: map_vector.hpp:211
Ghost
Definition: Ghost.hpp:39
petsc_solver< double >::setRestart
void setRestart(PetscInt n)
For GMRES based method, the number of Krylov directions to orthogonalize against.
Definition: petsc_solver.hpp:1014
Point::get
__device__ __host__ const T & get(unsigned int i) const
Get coordinate.
Definition: Point.hpp:172
petsc_solver< double >
This class is able to do Matrix inversion in parallel with PETSC solvers.
Definition: petsc_solver.hpp:101
petsc_solver< double >::setSolver
void setSolver(KSPType type)
Set the Petsc solver.
Definition: petsc_solver.hpp:918
Lap
Laplacian second order on h (spacing)
Definition: Laplacian.hpp:22
Box
This class represent an N-dimensional box.
Definition: Box.hpp:60
prop_id
Definition: eq_solve_common.hpp:12
vector_dist
Distributed vector.
Definition: vector_dist.hpp:297
openfpm::vector
Implementation of 1-D std::vector like structure.
Definition: map_vector.hpp:202
eq_id
Definition: eq_solve_common.hpp:14