OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
DCPSE_op_test3d.cpp
1 /*
2  * DCPSE_op_test.cpp
3  *
4  * Created on: April 9, 2020
5  * Author: Abhinav Singh
6  *
7  */
8 #include "config.h"
9 #ifdef HAVE_EIGEN
10 #ifdef HAVE_PETSC
11 
12 
13 
14 #define BOOST_TEST_DYN_LINK
15 
16 #include "util/util_debug.hpp"
17 #include <boost/test/unit_test.hpp>
18 #include <iostream>
19 #include "../DCPSE_op.hpp"
20 #include "../DCPSE_Solver.hpp"
21 #include "Operators/Vector/vector_dist_operators.hpp"
22 #include "Vector/vector_dist_subset.hpp"
23 #include "../EqnsStruct.hpp"
24 #include "util/SphericalHarmonics.hpp"
25 
26 
27 //template<typename T>
28 //struct Debug;
29 
30 
31 
32 BOOST_AUTO_TEST_SUITE(dcpse_op_suite_tests3)
33  BOOST_AUTO_TEST_CASE(dcpse_op_vec3d) {
34 // int rank;
35 // MPI_Comm_rank(MPI_COMM_WORLD, &rank);
36  size_t edgeSemiSize = 21;
37  const size_t sz[3] = {edgeSemiSize, edgeSemiSize,edgeSemiSize};
38  Box<3, double> box({0, 0,0}, {1,1,1});
39  size_t bc[3] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
40  double spacing = box.getHigh(0) / (sz[0] - 1);
41  double rCut = 3.1 * spacing;
42  Ghost<3, double> ghost(rCut);
43  BOOST_TEST_MESSAGE("Init vector_dist...");
44  double sigma2 = spacing * spacing/ (2 * 4);
45 
46  vector_dist<3, double, aggregate<double, VectorS<3, double>, VectorS<3, double>, VectorS<3, double>, VectorS<3, double>,double,double>> domain(
47  0, box, bc, ghost);
48 
49  //Init_DCPSE(domain)
50  BOOST_TEST_MESSAGE("Init domain...");
51 
52  auto it = domain.getGridIterator(sz);
53  size_t pointId = 0;
54  size_t counter = 0;
55  double minNormOne = 999;
56  while (it.isNext()) {
57  domain.add();
58  auto key = it.get();
59  mem_id k0 = key.get(0);
60  double x = k0 * spacing;
61  domain.getLastPos()[0] = x;//+ gaussian(rng);
62  mem_id k1 = key.get(1);
63  double y = k1 * spacing;
64  domain.getLastPos()[1] = y;//+gaussian(rng);
65  mem_id k2 = key.get(2);
66  double z = k2 * spacing;
67  domain.getLastPos()[2] = z;//+gaussian(rng);
68  // Here fill the function value
69  domain.template getLastProp<0>() = sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1]) + sin(domain.getLastPos()[2]) ;
70  domain.template getLastProp<1>()[0] = cos(domain.getLastPos()[0]);
71  domain.template getLastProp<1>()[1] = cos(domain.getLastPos()[1]) ;
72  domain.template getLastProp<1>()[2] = cos(domain.getLastPos()[2]);
73  // Here fill the validation value for Df/Dx
74  domain.template getLastProp<2>()[0] = 0;//cos(domain.getLastPos()[0]);//+cos(domain.getLastPos()[1]);
75  domain.template getLastProp<2>()[1] = 0;//-sin(domain.getLastPos()[0]);//+cos(domain.getLastPos()[1]);
76  domain.template getLastProp<3>()[0] = 0;//cos(domain.getLastPos()[0]);//+cos(domain.getLastPos()[1]);
77  domain.template getLastProp<3>()[1] = 0;//-sin(domain.getLastPos()[0]);//+cos(domain.getLastPos()[1]);
78  domain.template getLastProp<3>()[2] = 0;
79 
80  domain.template getLastProp<4>()[0] = -cos(domain.getLastPos()[0]) * sin(domain.getLastPos()[0]);
81  domain.template getLastProp<4>()[1] = -cos(domain.getLastPos()[1]) * sin(domain.getLastPos()[1]);
82  domain.template getLastProp<4>()[2] = -cos(domain.getLastPos()[2]) * sin(domain.getLastPos()[2]);
83 
84 
85  /* domain.template getLastProp<4>()[0] = cos(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) +
86  cos(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
87  domain.template getLastProp<4>()[1] = -sin(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) -
88  sin(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));
89  domain.template getLastProp<4>()[2] = -sin(domain.getLastPos()[0]) * (sin(domain.getLastPos()[0]) + sin(domain.getLastPos()[1])) -
90  sin(domain.getLastPos()[1]) * (cos(domain.getLastPos()[0]) + cos(domain.getLastPos()[1]));*/
91  domain.template getLastProp<5>() = cos(domain.getLastPos()[0]) * cos(domain.getLastPos()[0])+cos(domain.getLastPos()[1]) * cos(domain.getLastPos()[1])+cos(domain.getLastPos()[2]) * cos(domain.getLastPos()[2]) ;
92  ++counter;
93  ++it;
94  }
95  BOOST_TEST_MESSAGE("Sync domain across processors...");
96 
97  domain.map();
98  domain.ghost_get<0>();
99 
100  Advection Adv(domain, 2, rCut, 1.9,support_options::RADIUS);
101  auto v = getV<1>(domain);
102  auto P = getV<0>(domain);
103  auto dv = getV<3>(domain);
104  auto dP = getV<6>(domain);
105 
106 
107 // typedef boost::mpl::int_<std::is_fundamental<point_expression_op<Point<2U, double>, point_expression<double>, Point<2U, double>, 3>>::value>::blabla blabla;
108 // std::is_fundamental<decltype(o1.value(key))>
109 
110  domain.ghost_get<1>();
111  dv = Adv(v, v);
112  auto it2 = domain.getDomainIterator();
113 
114  double worst1 = 0.0;
115 
116  while (it2.isNext()) {
117  auto p = it2.get();
118 
119  if (fabs(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]) > worst1) {
120  worst1 = fabs(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]);
121 
122  }
123 
124  ++it2;
125  }
126  //std::cout << "Maximum Error in component 2: " << worst1 << std::endl;
127  BOOST_REQUIRE(worst1 < 0.03);
128 
129  //Adv.checkMomenta(domain);
130  //Adv.DrawKernel<2>(domain,0);
131 
132  //domain.deleteGhost();
133 
134  dP = Adv(v, P);//+Dy(P);
135  auto it3 = domain.getDomainIterator();
136 
137  double worst2 = 0.0;
138 
139  while (it3.isNext()) {
140  auto p = it3.get();
141  if (fabs(domain.getProp<6>(p) - domain.getProp<5>(p)) > worst2) {
142  worst2 = fabs(domain.getProp<6>(p) - domain.getProp<5>(p));
143 
144  }
145 
146  ++it3;
147  }
148  domain.deleteGhost();
149  BOOST_REQUIRE(worst2 < 0.03);
150 
151 
152  }
154  BOOST_AUTO_TEST_CASE(dcpse_poisson_dirichlet_anal3d) {
155 // int rank;
156 // MPI_Comm_rank(MPI_COMM_WORLD, &rank);
157  size_t grd_sz=21;
158  const size_t sz[3] = {grd_sz,grd_sz,grd_sz};
159  Box<3, double> box({0, 0,0}, {1.0, 1.0,1.0});
160  size_t bc[3] = {NON_PERIODIC, NON_PERIODIC,NON_PERIODIC};
161  double spacing = box.getHigh(0) / (sz[0] - 1);
162  Ghost<3, double> ghost(spacing * 3.1);
163  double rCut = 3.1 * spacing;
164  BOOST_TEST_MESSAGE("Init vector_dist...");
165 
166  vector_dist<3, double, aggregate<double,double,double,double,double,double>> domain(0, box, bc, ghost);
167 
168 
169  //Init_DCPSE(domain)
170  BOOST_TEST_MESSAGE("Init domain...");
171 
172  auto it = domain.getGridIterator(sz);
173  while (it.isNext()) {
174  domain.add();
175  auto key = it.get();
176  double x = key.get(0) * it.getSpacing(0);
177  domain.getLastPos()[0] = x;
178  double y = key.get(1) * it.getSpacing(1);
179  domain.getLastPos()[1] = y;
180  double z = key.get(2) * it.getSpacing(2);
181  domain.getLastPos()[2] = z;
182 
183  ++it;
184  }
185  BOOST_TEST_MESSAGE("Sync domain across processors...");
186 
187  domain.map();
188  domain.ghost_get<0>();
189 
190  Derivative_x Dx(domain, 2, rCut,1.9,support_options::RADIUS);
191  Derivative_y Dy(domain, 2, rCut,1.9,support_options::RADIUS);
192  Laplacian Lap(domain, 2, rCut,1.3,support_options::RADIUS);
193 
194  openfpm::vector<aggregate<int>> bulk;
195  openfpm::vector<aggregate<int>> front_p;
196  openfpm::vector<aggregate<int>> back_p;
197  openfpm::vector<aggregate<int>> up_p;
198  openfpm::vector<aggregate<int>> down_p;
199  openfpm::vector<aggregate<int>> left_p;
200  openfpm::vector<aggregate<int>> right_p;
201 
202  auto v = getV<0>(domain);
203  auto RHS=getV<1>(domain);
204  auto sol = getV<2>(domain);
205  auto anasol = getV<3>(domain);
206  auto err = getV<4>(domain);
207  auto DCPSE_sol=getV<5>(domain);
208 
209  Box<3, double> up(
210  {box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0, box.getLow(2) - spacing / 2.0},
211  {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0, box.getHigh(2) + spacing / 2.0});
212 
213  Box<3, double> down(
214  {box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0, box.getLow(2) - spacing / 2.0},
215  {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0, box.getHigh(2) + spacing / 2.0});
216 
217  Box<3, double> left(
218  {box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0, box.getLow(2) - spacing / 2.0},
219  {box.getLow(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0, box.getHigh(2) + spacing / 2.0});
220 
221  Box<3, double> right(
222  {box.getHigh(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0, box.getLow(2) - spacing / 2.0},
223  {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0, box.getHigh(2) + spacing / 2.0});
224 
225  Box<3, double> front(
226  {box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0, box.getLow(2) - spacing / 2.0},
227  {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0, box.getLow(2) + spacing / 2.0});
228 
229  Box<3, double> back(
230  {box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0, box.getHigh(2) - spacing / 2.0},
231  {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0, box.getHigh(2) + spacing / 2.0});
232 
233  openfpm::vector<Box<3, double>> boxes;
234  boxes.add(up);
235  boxes.add(down);
236  boxes.add(left);
237  boxes.add(right);
238  boxes.add(front);
239  boxes.add(back);
240  VTKWriter<openfpm::vector<Box<3, double>>, VECTOR_BOX> vtk_box;
241  vtk_box.add(boxes);
242  //vtk_box.write("boxes_3d.vtk");
243  auto Particles=domain;
244  auto it2 = Particles.getDomainIterator();
245 
246  while (it2.isNext()) {
247  auto p = it2.get();
248  Point<3, double> xp = Particles.getPos(p);
249  domain.getProp<1>(p) = -3.0*M_PI*M_PI*sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1))*sin(M_PI*xp.get(2));
250  domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1))*sin(M_PI*xp.get(2));
251  if (front.isInside(xp) == true) {
252  front_p.add();
253  front_p.last().get<0>() = p.getKey();
254  } else if (back.isInside(xp) == true) {
255  back_p.add();
256  back_p.last().get<0>() = p.getKey();
257  } else if (left.isInside(xp) == true) {
258  left_p.add();
259  left_p.last().get<0>() = p.getKey();
260  } else if (right.isInside(xp) == true) {
261  right_p.add();
262  right_p.last().get<0>() = p.getKey();
263  } else if (up.isInside(xp) == true) {
264  up_p.add();
265  up_p.last().get<0>() = p.getKey();
266  } else if (down.isInside(xp) == true) {
267  down_p.add();
268  down_p.last().get<0>() = p.getKey();
269  } else {
270  bulk.add();
271  bulk.last().get<0>() = p.getKey();
272  }
273  ++it2;
274  }
275 
276 
277  DCPSE_scheme<equations3d1,decltype(domain)> Solver( domain);
278  auto Poisson = Lap(v);
279  auto D_x = Dx(v);
280  auto D_y = Dy(v);
281  Solver.impose(Poisson, bulk, prop_id<1>());
282  Solver.impose(v, up_p, prop_id<1>());
283  Solver.impose(v, right_p, prop_id<1>());
284  Solver.impose(v, down_p, prop_id<1>());
285  Solver.impose(v, left_p, prop_id<1>());
286  Solver.impose(v, front_p, prop_id<1>());
287  Solver.impose(v, back_p, prop_id<1>());
288  Solver.solve(sol);
289  DCPSE_sol=Lap(sol);
290 
291  double worst1 = 0.0;
292 
293  v=abs(DCPSE_sol-RHS);
294 
295  for(int j=0;j<bulk.size();j++)
296  { auto p=bulk.get<0>(j);
297  if (fabs(domain.getProp<3>(p) - domain.getProp<2>(p)) >= worst1) {
298  worst1 = fabs(domain.getProp<3>(p) - domain.getProp<2>(p));
299  }
300  domain.getProp<4>(p) = fabs(domain.getProp<3>(p) - domain.getProp<2>(p));
301 
302  }
303  //std::cout << "Maximum Analytic Error: " << worst1 << std::endl;
304  //domain.write("Dirichlet_anasol_3d");
305 
306  BOOST_REQUIRE(worst1 < 0.031);
307  }
308 
309 //Is failing on Ubuntu CI with 5 cores. Needs investigation.
310 /* BOOST_AUTO_TEST_CASE(Sph_harm) {
311  BOOST_REQUIRE(openfpm::math::Y(2,1,0.5,0)+0.459674<0.00001);
312  //These would be a requirement once Boost releases their fix
313  //
314  //BOOST_REQUIRE(boost::math::legendre_p(0,-1,1)=?);
315  double nu=1.0;
316  size_t grd_sz=20;
317  const size_t sz[3] = {grd_sz,grd_sz,grd_sz};
318  Box<3, double> box({-1.0, -1.0,-1.0}, {1.0,1.0,1.0});
319  size_t bc[3] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
320  double spacing = 2.0 / (sz[0] - 1);
321  double rCut = 3.9*spacing;
322  double R=1.0;
323  Ghost<3, double> ghost(rCut);
324  // P V v_B RHS V_t P_anal RHS2 Polar cord
325  vector_dist_ws<3, double, aggregate<double,VectorS<3, double>,VectorS<3, double>,double,VectorS<3, double>,double,double,VectorS<3, double>,VectorS<3, double>,VectorS<3, double>>> Particles(0, box, bc, ghost);
326 
327 
328  auto &v_cl = create_vcluster();
329 
330 // openfpm::vector<aggregate<int>> bulk;
331 // openfpm::vector<aggregate<int>> Surface;
332 
333  auto it = Particles.getGridIterator(sz);
334  while (it.isNext()) {
335  auto key = it.get();
336  double x = -1.0+key.get(0) * it.getSpacing(0);
337  double y = -1.0+key.get(1) * it.getSpacing(1);
338  double z = -1.0+key.get(2) * it.getSpacing(2);
339  double r=sqrt(x*x+y*y+z*z);
340  if (r<R-spacing/2.0) {
341  Particles.add();
342  Particles.getLastPos()[0] = x;
343  Particles.getLastPos()[1] = y;
344  Particles.getLastPos()[2] = z;
345  Particles.getLastProp<8>()[0] = r;
346  if (r==0){
347  Particles.getLastProp<8>()[1] = 0.0;
348  }
349  else{
350  Particles.getLastProp<8>()[1] = std::atan2(sqrt(x*x+y*y),z);
351  }
352  Particles.getLastProp<8>()[2] = std::atan2(y,x);
353  }
354  ++it;
355  }
356 
357  int n_sp=int(grd_sz)*int(grd_sz)*3;
358 
359  double Golden_angle=M_PI * (3.0 - sqrt(5.0));
360 
361  for(int i=1;i<=n_sp;i++)
362  {
363  double y = 1.0 - (i /double(n_sp - 1.0)) * 2.0;
364  double radius = sqrt(1 - y * y);
365  double Golden_theta = Golden_angle * i;
366  double x = cos(Golden_theta) * radius;
367  double z = sin(Golden_theta) * radius;
368 
369  if (acos(z)==0 || acos(z)==M_PI){
370  std::cout<<"Theta 0/Pi "<<std::endl;
371  continue;
372  }
373 
374  Particles.add();
375  Particles.getLastPos()[0] = x;
376  Particles.getLastPos()[1] = y;
377  Particles.getLastPos()[2] = z;
378  Particles.getLastProp<8>()[0] = 1.0 ;
379  Particles.getLastProp<8>()[1] = std::atan2(sqrt(x*x+y*y),z);
380  Particles.getLastProp<8>()[2] = std::atan2(y,x);
381  }
382  Particles.map();
383  Particles.ghost_get<0>();
384 
385 
386  std::unordered_map<const lm,double,key_hash,key_equal> Vr;
387  std::unordered_map<const lm,double,key_hash,key_equal> V1;
388  std::unordered_map<const lm,double,key_hash,key_equal> V2;
389  //Setting max mode l_max
390  constexpr int K = 2;
391  //Setting amplitudes to 0
392  for(int l=0;l<=K;l++){
393  for(int m=-l;m<=l;m++){
394  Vr[std::make_tuple(l,m)]=0.0;
395  V1[std::make_tuple(l,m)]=0.0;
396  V2[std::make_tuple(l,m)]=0.0;
397  }
398 
399 
400  }
401  //Setting some amplitude for boundary velocity
402  V1[std::make_tuple(1,0)]=1.0;
403 
404  auto it2 = Particles.getDomainIterator();
405  while (it2.isNext()) {
406  auto p = it2.get();
407  Point<3, double> xp = Particles.getPos(p);
408  Point<3, double> xP = Particles.getProp<8>(p);
409  Particles.getProp<0>(p) =0;
410  if (xP[0]==1.0) {
411 // Surface.add();
412 // Surface.last().get<0>() = p.getKey();
413  Particles.getProp<0>(p) = 0;
414  std::vector<double> SVel;
415  SVel=openfpm::math::sumY<K>(xP[0],xP[1],xP[2],Vr,V1,V2);
416  double SP=openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Vr);
417  Particles.getProp<2>(p)[0] = SVel[0];
418  Particles.getProp<2>(p)[1] = SVel[1];
419  Particles.getProp<2>(p)[2] = SVel[2];
420  Particles.getProp<9>(p)[0] = SVel[0];
421  Particles.getProp<9>(p)[1] = SVel[1];
422  Particles.getProp<9>(p)[2] = SVel[2];
423  Particles.getProp<5>(p) = SP;
424  Particles.setSubset(p,1);
425 
426  }
427  else {
428 // bulk.add();
429 // bulk.last().get<0>() = p.getKey();
430  Particles.setSubset(p,0);
431  Particles.getProp<0>(p) = 0;
432  Particles.getProp<1>(p)[0] = 0;
433  Particles.getProp<1>(p)[1] = 0;
434  Particles.getProp<1>(p)[2] = 0;
435  }
436  ++it2;
437  }
438 
439  vector_dist_subset<3, double, aggregate<double,VectorS<3, double>,VectorS<3, double>,double,VectorS<3, double>,double,double,VectorS<3, double>,VectorS<3, double>,VectorS<3, double>>> Particles_bulk(Particles,0);
440  vector_dist_subset<3, double, aggregate<double,VectorS<3, double>,VectorS<3, double>,double,VectorS<3, double>,double,double,VectorS<3, double>,VectorS<3, double>,VectorS<3, double>>> Particles_surface(Particles,1);
441 
442  auto & bulk = Particles_bulk.getIds();
443  auto & Surface = Particles_surface.getIds();
444 
445  for (int j = 0; j < bulk.size(); j++) {
446  auto p = bulk.get<0>(j);
447  Point<3, double> xp = Particles.getPos(p);
448  Point<3, double> xP = Particles.getProp<8>(p);
449 
450  std::unordered_map<const lm,double,key_hash,key_equal> Ur;
451  std::unordered_map<const lm,double,key_hash,key_equal> U2;
452  std::unordered_map<const lm,double,key_hash,key_equal> U1;
453  std::unordered_map<const lm,double,key_hash,key_equal> Plm;
454 
455  for (int l = 0; l <= K; l++) {
456  for (int m = -l; m <= l; m++) {
457  auto Er= Vr.find(std::make_tuple(l,m));
458  auto E1= V1.find(std::make_tuple(l,m));
459  auto E2= V2.find(std::make_tuple(l,m));
460  std::vector<double> Sol=openfpm::math::sph_anasol_u(nu,l,m,Er->second,E1->second,E2->second,xP[0]);
461  Ur[std::make_tuple(l,m)]=Sol[0];
462  U1[std::make_tuple(l,m)]=Sol[1];
463  U2[std::make_tuple(l,m)]=Sol[2];
464  Plm[std::make_tuple(l,m)]=Sol[3];
465  }
466 
467  }
468 
469  if(fabs(xP[0])>=1e-5 && xP[1]>1e-5 && (M_PI-xP[1])>=1e-5)
470  {
471  std::vector<double> SVel = openfpm::math::sumY<K>(xP[0], xP[1], xP[2], Ur, U1, U2);
472  Particles.getProp<9>(p)[0] = SVel[0];
473  Particles.getProp<9>(p)[1] = SVel[1];
474  Particles.getProp<9>(p)[2] = SVel[2];
475  Particles.getProp<5>(p) = openfpm::math::sumY_Scalar<K>(xP[0], xP[1], xP[2], Plm);
476  }
477  }
478 
479  auto P = getV<0>(Particles);
480  auto V = getV<1>(Particles);
481  auto V_B = getV<2>(Particles);
482  V.setVarId(0);
483  auto DIV = getV<3>(Particles);
484  auto V_t = getV<4>(Particles);
485  auto P_anal = getV<5>(Particles);
486  auto temp=getV<6>(Particles);
487  auto RHS = getV<7>(Particles);
488  auto P_bulk = getV<0>(Particles_bulk);
489  auto RHS_bulk = getV<7>(Particles_bulk);
490  auto V_anal = getV<9>(Particles);
491 
492  V_t=V;
493  P=0;
494  P_bulk=0;
495  eq_id vx,vy,vz;
496 
497  vx.setId(0);
498  vy.setId(1);
499  vz.setId(2);
500 
501  double sampling=3.1;
502  double sampling2=1.9;
503  double rCut2=3.9*spacing;
504 
505  Derivative_x Dx(Particles, 2, rCut,sampling, support_options::RADIUS),B_Dx(Particles_bulk, 2, rCut,sampling, support_options::RADIUS);
506  Derivative_y Dy(Particles, 2, rCut,sampling, support_options::RADIUS),B_Dy(Particles_bulk, 2, rCut,sampling, support_options::RADIUS);
507  Derivative_z Dz(Particles, 2, rCut,sampling, support_options::RADIUS),B_Dz(Particles_bulk, 2, rCut,sampling, support_options::RADIUS);
508  Derivative_xx Dxx(Particles, 2, rCut2,sampling2,support_options::RADIUS);
509  Derivative_yy Dyy(Particles, 2, rCut2,sampling2,support_options::RADIUS);
510  Derivative_zz Dzz(Particles, 2, rCut2,sampling2,support_options::RADIUS);
511 
512  //std::cout << "DCPSE KERNELS DONE" << std::endl;
513  petsc_solver<double> solverPetsc;
514  solverPetsc.setPreconditioner(PCNONE);
515  timer tt;
516  double sum=0,sum1=0;
517  V_t=V;
518  double V_err_eps = 1e-5;
519 
520  double V_err = 1, V_err_old;
521  int n = 0;
522  int nmax = 30;
523  int ctr = 0, errctr, Vreset = 0;
524  V_err = 1;
525  n = 0;
526  tt.start();
527  while (V_err >= V_err_eps && n <= nmax) {
528  //Particles.write_frame("StokesSphere",n);
529  Particles.ghost_get<0>(SKIP_LABELLING);
530  RHS_bulk[0] = B_Dx(P);
531  RHS_bulk[1] = B_Dy(P);
532  RHS_bulk[2] = B_Dz(P);
533  DCPSE_scheme<equations3d3, decltype(Particles)> Solver(Particles);
534  auto Stokes1 = nu * (Dxx(V[0])+Dyy(V[0])+Dzz(V[0]));
535  auto Stokes2 = nu * (Dxx(V[1])+Dyy(V[1])+Dzz(V[1]));
536  auto Stokes3 = nu * (Dxx(V[2])+Dyy(V[2])+Dzz(V[2]));
537  Solver.impose(Stokes1, bulk, RHS[0], vx);
538  Solver.impose(Stokes2, bulk, RHS[1], vy);
539  Solver.impose(Stokes3, bulk, RHS[2], vz);
540  Solver.impose(V[0], Surface, V_B[0], vx);
541  Solver.impose(V[1], Surface, V_B[1], vy);
542  Solver.impose(V[2], Surface, V_B[2], vz);
543  Solver.solve_with_solver(solverPetsc, V[0], V[1], V[2]);
544  //Solver.solve(V[0],V[1],V[2]);
545  //std::cout << "Stokes Solved" << std::endl;
546  Particles.ghost_get<1>();
547  DIV = -(Dx(V[0])+Dy(V[1])+Dz(V[2]));
548  P_bulk = P + DIV;
549  sum = 0;
550  sum1 = 0;
551  for (int j = 0; j < bulk.size(); j++) {
552  auto p = bulk.get<0>(j);
553  sum += (Particles.getProp<4>(p)[0] - Particles.getProp<1>(p)[0]) *
554  (Particles.getProp<4>(p)[0] - Particles.getProp<1>(p)[0]) +
555  (Particles.getProp<4>(p)[1] - Particles.getProp<1>(p)[1]) *
556  (Particles.getProp<4>(p)[1] - Particles.getProp<1>(p)[1]) +
557  (Particles.getProp<4>(p)[2] - Particles.getProp<1>(p)[2]) *
558  (Particles.getProp<4>(p)[2] - Particles.getProp<1>(p)[2]);
559  sum1 += Particles.getProp<1>(p)[0] * Particles.getProp<1>(p)[0] +
560  Particles.getProp<1>(p)[1] * Particles.getProp<1>(p)[1] +
561  Particles.getProp<1>(p)[2] * Particles.getProp<1>(p)[2];
562  }
563  sum = sqrt(sum);
564  sum1 = sqrt(sum1);
565  v_cl.sum(sum);
566  v_cl.sum(sum1);
567  v_cl.execute();
568  V_t = V;
569  Particles.ghost_get<1>(SKIP_LABELLING);
570  V_err_old = V_err;
571  V_err = sum / sum1;
572  if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-14) {
573  errctr++;
574  } else {
575  errctr = 0;
576  }
577  if (n > 3) {
578  if (errctr > 1) {
579  std::cout << "CONVERGENCE LOOP BROKEN DUE TO INCREASE/VERY SLOW DECREASE IN ERROR" << std::endl;
580  Vreset = 1;
581  break;
582  } else {
583  Vreset = 0;
584  }
585  }
586  n++;
587 
588  }
589 
590  V_t=0;
591 
592  double worst=0;
593  double L2=0;
594  for (int j = 0; j < bulk.size(); j++) {
595  auto p = bulk.get<0>(j);
596  Point<3,double> xP=Particles.getProp<8>(p);
597  if(xP[0]>=1e-5 && xP[1]>1e-5 && (M_PI-xP[1])>=1e-5)
598  {
599  double dx=Particles.getProp<9>(p)[0] - Particles.getProp<1>(p)[0];
600  double dy=Particles.getProp<9>(p)[1] - Particles.getProp<1>(p)[1];
601  double dz=Particles.getProp<9>(p)[2] - Particles.getProp<1>(p)[2];
602  Particles.getProp<4>(p)[0]=fabs(dx);
603  Particles.getProp<4>(p)[1]=fabs(dy);
604  Particles.getProp<4>(p)[2]=fabs(dz);
605  L2 += dx*dx+dy*dy+dz*dz;
606  if (std::max({fabs(dx),fabs(dy),fabs(dz)}) > worst) {
607  worst = std::max({fabs(dx),fabs(dy),fabs(dz)});
608  }
609 
610  }
611  }
612 
613  v_cl.sum(worst);
614  v_cl.sum(L2);
615  v_cl.execute();
616 *//* if (v_cl.rank() == 0) {
617  std::cout<<"Gd,Surf,Bulk Size: "<<grd_sz<<","<<Surface.size()<<","<<bulk.size()<<std::endl;
618  std::cout << "L2_Final: " <<sqrt(L2)<<","<<sqrt(L2/(bulk.size()+Surface.size()))
619  << std::endl;
620  std::cout << "L_inf_Final: " << worst
621  << std::endl;
622  }*//*
623  std::cout << "L_inf_Final_test: " << worst;
624  //Particles.write("StokesSphere");
625  BOOST_REQUIRE(worst<1e-3);
626 
627  }*/
628 
629 BOOST_AUTO_TEST_SUITE_END()
630 
631 #endif
632 #endif
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
Point::get
__device__ __host__ const T & get(unsigned int i) const
Get coordinate.
Definition: Point.hpp:172
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
Point_test
Test structure used for several test.
Definition: Point_test.hpp:105
openfpm::vector
Implementation of 1-D std::vector like structure.
Definition: map_vector.hpp:202