OpenFPM  5.2.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/DCPSE_op/DCPSE_op.hpp"
20 #include "DCPSE/DCPSE_op/DCPSE_Solver.hpp"
21 #include "Operators/Vector/vector_dist_operators.hpp"
22 #include "Vector/vector_dist_subset.hpp"
23 #include "DCPSE/DCPSE_op/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  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
101  Advection<decltype(verletList)> Adv(domain, verletList, 2, rCut, support_options::RADIUS);
102  auto v = getV<1>(domain);
103  auto P = getV<0>(domain);
104  auto dv = getV<3>(domain);
105  auto dP = getV<6>(domain);
106 
107 
108 // typedef boost::mpl::int_<std::is_fundamental<point_expression_op<Point<2U, double>, point_expression<double>, Point<2U, double>, 3>>::value>::blabla blabla;
109 // std::is_fundamental<decltype(o1.value(key))>
110 
111  domain.ghost_get<1>();
112  dv = Adv(v, v);
113  auto it2 = domain.getDomainIterator();
114 
115  double worst1 = 0.0;
116 
117  while (it2.isNext()) {
118  auto p = it2.get();
119 
120  if (fabs(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]) > worst1) {
121  worst1 = fabs(domain.getProp<3>(p)[1] - domain.getProp<4>(p)[1]);
122 
123  }
124 
125  ++it2;
126  }
127  //std::cout << "Maximum Error in component 2: " << worst1 << std::endl;
128  BOOST_REQUIRE(worst1 < 0.03);
129 
130  //Adv.checkMomenta(domain);
131  //Adv.DrawKernel<2>(domain,0);
132 
133  //domain.deleteGhost();
134 
135  dP = Adv(v, P);//+Dy(P);
136  auto it3 = domain.getDomainIterator();
137 
138  double worst2 = 0.0;
139 
140  while (it3.isNext()) {
141  auto p = it3.get();
142  if (fabs(domain.getProp<6>(p) - domain.getProp<5>(p)) > worst2) {
143  worst2 = fabs(domain.getProp<6>(p) - domain.getProp<5>(p));
144 
145  }
146 
147  ++it3;
148  }
149  domain.deleteGhost();
150  BOOST_REQUIRE(worst2 < 0.03);
151 
152 
153  }
155  BOOST_AUTO_TEST_CASE(dcpse_poisson_dirichlet_anal3d) {
156 // int rank;
157 // MPI_Comm_rank(MPI_COMM_WORLD, &rank);
158  size_t grd_sz=21;
159  const size_t sz[3] = {grd_sz,grd_sz,grd_sz};
160  Box<3, double> box({0, 0,0}, {1.0, 1.0,1.0});
161  size_t bc[3] = {NON_PERIODIC, NON_PERIODIC,NON_PERIODIC};
162  double spacing = box.getHigh(0) / (sz[0] - 1);
163  Ghost<3, double> ghost(spacing * 3.1);
164  double rCut = 3.1 * spacing;
165  BOOST_TEST_MESSAGE("Init vector_dist...");
166 
167  vector_dist<3, double, aggregate<double,double,double,double,double,double>> domain(0, box, bc, ghost);
168 
169 
170  //Init_DCPSE(domain)
171  BOOST_TEST_MESSAGE("Init domain...");
172 
173  auto it = domain.getGridIterator(sz);
174  while (it.isNext()) {
175  domain.add();
176  auto key = it.get();
177  double x = key.get(0) * it.getSpacing(0);
178  domain.getLastPos()[0] = x;
179  double y = key.get(1) * it.getSpacing(1);
180  domain.getLastPos()[1] = y;
181  double z = key.get(2) * it.getSpacing(2);
182  domain.getLastPos()[2] = z;
183 
184  ++it;
185  }
186  BOOST_TEST_MESSAGE("Sync domain across processors...");
187 
188  domain.map();
189  domain.ghost_get<0>();
190 
191  auto verletList = domain.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
192 
193  Derivative_x<decltype(verletList)> Dx(domain, verletList, 2, rCut, support_options::RADIUS);
194  Derivative_y<decltype(verletList)> Dy(domain, verletList, 2, rCut, support_options::RADIUS);
195  Laplacian<decltype(verletList)> Lap(domain, verletList, 2, rCut, support_options::RADIUS);
196 
197  openfpm::vector<aggregate<int>> bulk;
198  openfpm::vector<aggregate<int>> front_p;
199  openfpm::vector<aggregate<int>> back_p;
200  openfpm::vector<aggregate<int>> up_p;
201  openfpm::vector<aggregate<int>> down_p;
202  openfpm::vector<aggregate<int>> left_p;
203  openfpm::vector<aggregate<int>> right_p;
204 
205  auto v = getV<0>(domain);
206  auto RHS=getV<1>(domain);
207  auto sol = getV<2>(domain);
208  auto anasol = getV<3>(domain);
209  auto err = getV<4>(domain);
210  auto DCPSE_sol=getV<5>(domain);
211 
212  Box<3, double> up(
213  {box.getLow(0) - spacing / 2.0, box.getHigh(1) - spacing / 2.0, box.getLow(2) - spacing / 2.0},
214  {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0, box.getHigh(2) + spacing / 2.0});
215 
216  Box<3, double> down(
217  {box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0, box.getLow(2) - spacing / 2.0},
218  {box.getHigh(0) + spacing / 2.0, box.getLow(1) + spacing / 2.0, box.getHigh(2) + spacing / 2.0});
219 
220  Box<3, double> left(
221  {box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0, box.getLow(2) - spacing / 2.0},
222  {box.getLow(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0, box.getHigh(2) + spacing / 2.0});
223 
224  Box<3, double> right(
225  {box.getHigh(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0, box.getLow(2) - spacing / 2.0},
226  {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0, box.getHigh(2) + spacing / 2.0});
227 
228  Box<3, double> front(
229  {box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0, box.getLow(2) - spacing / 2.0},
230  {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0, box.getLow(2) + spacing / 2.0});
231 
232  Box<3, double> back(
233  {box.getLow(0) - spacing / 2.0, box.getLow(1) - spacing / 2.0, box.getHigh(2) - spacing / 2.0},
234  {box.getHigh(0) + spacing / 2.0, box.getHigh(1) + spacing / 2.0, box.getHigh(2) + spacing / 2.0});
235 
236  openfpm::vector<Box<3, double>> boxes;
237  boxes.add(up);
238  boxes.add(down);
239  boxes.add(left);
240  boxes.add(right);
241  boxes.add(front);
242  boxes.add(back);
243  VTKWriter<openfpm::vector<Box<3, double>>, VECTOR_BOX> vtk_box;
244  vtk_box.add(boxes);
245  //vtk_box.write("boxes_3d.vtk");
246  auto Particles=domain;
247  auto it2 = Particles.getDomainIterator();
248 
249  while (it2.isNext()) {
250  auto p = it2.get();
251  Point<3, double> xp = Particles.getPos(p);
252  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));
253  domain.getProp<3>(p) = sin(M_PI*xp.get(0))*sin(M_PI*xp.get(1))*sin(M_PI*xp.get(2));
254  if (front.isInside(xp) == true) {
255  front_p.add();
256  front_p.last().get<0>() = p.getKey();
257  } else if (back.isInside(xp) == true) {
258  back_p.add();
259  back_p.last().get<0>() = p.getKey();
260  } else if (left.isInside(xp) == true) {
261  left_p.add();
262  left_p.last().get<0>() = p.getKey();
263  } else if (right.isInside(xp) == true) {
264  right_p.add();
265  right_p.last().get<0>() = p.getKey();
266  } else if (up.isInside(xp) == true) {
267  up_p.add();
268  up_p.last().get<0>() = p.getKey();
269  } else if (down.isInside(xp) == true) {
270  down_p.add();
271  down_p.last().get<0>() = p.getKey();
272  } else {
273  bulk.add();
274  bulk.last().get<0>() = p.getKey();
275  }
276  ++it2;
277  }
278 
279 
280  DCPSE_scheme<equations3d1,decltype(domain)> Solver( domain);
281  auto Poisson = Lap(v);
282  auto D_x = Dx(v);
283  auto D_y = Dy(v);
284  Solver.impose(Poisson, bulk, prop_id<1>());
285  Solver.impose(v, up_p, prop_id<1>());
286  Solver.impose(v, right_p, prop_id<1>());
287  Solver.impose(v, down_p, prop_id<1>());
288  Solver.impose(v, left_p, prop_id<1>());
289  Solver.impose(v, front_p, prop_id<1>());
290  Solver.impose(v, back_p, prop_id<1>());
291  Solver.solve(sol);
292  DCPSE_sol=Lap(sol);
293 
294  double worst1 = 0.0;
295 
296  v=abs(DCPSE_sol-RHS);
297 
298  for(int j=0;j<bulk.size();j++)
299  { auto p=bulk.get<0>(j);
300  if (fabs(domain.getProp<3>(p) - domain.getProp<2>(p)) >= worst1) {
301  worst1 = fabs(domain.getProp<3>(p) - domain.getProp<2>(p));
302  }
303  domain.getProp<4>(p) = fabs(domain.getProp<3>(p) - domain.getProp<2>(p));
304 
305  }
306  //std::cout << "Maximum Analytic Error: " << worst1 << std::endl;
307  //domain.write("Dirichlet_anasol_3d");
308 
309  BOOST_REQUIRE(worst1 < 0.031);
310  }
311 
312  BOOST_AUTO_TEST_CASE(Sph_harm) {
313  BOOST_REQUIRE(openfpm::math::Y(2,1,0.5,0)+0.459674<0.00001);
314  //These would be a requirement once Boost releases their fix
315  //
316  //BOOST_REQUIRE(boost::math::legendre_p(0,-1,1)=?);
317  double nu=1.0;
318  size_t grd_sz=13;
319  const size_t sz[3] = {grd_sz,grd_sz,grd_sz};
320  Box<3, double> box({-1.0, -1.0,-1.0}, {1.0,1.0,1.0});
321  size_t bc[3] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
322  double spacing = 2.0 / (sz[0] - 1);
323  double rCut = 3.9*spacing;
324  double R=1.0;
325  Ghost<3, double> ghost(rCut);
326  // P V v_B RHS V_t P_anal RHS2 Polar cord
327  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);
328 
329 
330  auto &v_cl = create_vcluster();
331 
332 // openfpm::vector<aggregate<int>> bulk;
333 // openfpm::vector<aggregate<int>> Surface;
334 
335  auto it = Particles.getGridIterator(sz);
336  while (it.isNext()) {
337  auto key = it.get();
338  double x = -1.0+key.get(0) * it.getSpacing(0);
339  double y = -1.0+key.get(1) * it.getSpacing(1);
340  double z = -1.0+key.get(2) * it.getSpacing(2);
341  double r=sqrt(x*x+y*y+z*z);
342  if (r<R-spacing/2.0) {
343  Particles.add();
344  Particles.getLastPos()[0] = x;
345  Particles.getLastPos()[1] = y;
346  Particles.getLastPos()[2] = z;
347  Particles.getLastProp<8>()[0] = r;
348  if (r==0){
349  Particles.getLastProp<8>()[1] = 0.0;
350  }
351  else{
352  Particles.getLastProp<8>()[1] = std::atan2(sqrt(x*x+y*y),z);
353  }
354  Particles.getLastProp<8>()[2] = std::atan2(y,x);
355  }
356  ++it;
357  }
358 
359  int n_sp=int(grd_sz)*int(grd_sz)*3;
360 
361  double Golden_angle=M_PI * (3.0 - sqrt(5.0));
362 
363  for(int i=1;i<=n_sp;i++)
364  {
365  double y = 1.0 - (i /double(n_sp - 1.0)) * 2.0;
366  double radius = sqrt(1 - y * y);
367  double Golden_theta = Golden_angle * i;
368  double x = cos(Golden_theta) * radius;
369  double z = sin(Golden_theta) * radius;
370 
371  if (acos(z)==0 || acos(z)==M_PI){
372  std::cout<<"Theta 0/Pi "<<std::endl;
373  continue;
374  }
375 
376  Particles.add();
377  Particles.getLastPos()[0] = x;
378  Particles.getLastPos()[1] = y;
379  Particles.getLastPos()[2] = z;
380  Particles.getLastProp<8>()[0] = 1.0 ;
381  Particles.getLastProp<8>()[1] = std::atan2(sqrt(x*x+y*y),z);
382  Particles.getLastProp<8>()[2] = std::atan2(y,x);
383  }
384  Particles.map();
385  Particles.ghost_get<0>();
386 
387 
388  std::unordered_map<const lm,double,key_hash,key_equal> Vr;
389  std::unordered_map<const lm,double,key_hash,key_equal> V1;
390  std::unordered_map<const lm,double,key_hash,key_equal> V2;
391  //Setting max mode l_max
392  constexpr int K = 2;
393  //Setting amplitudes to 0
394  for(int l=0;l<=K;l++){
395  for(int m=-l;m<=l;m++){
396  Vr[std::make_tuple(l,m)]=0.0;
397  V1[std::make_tuple(l,m)]=0.0;
398  V2[std::make_tuple(l,m)]=0.0;
399  }
400 
401 
402  }
403  //Setting some amplitude for boundary velocity
404  V1[std::make_tuple(1,0)]=1.0;
405 
406  auto it2 = Particles.getDomainIterator();
407  while (it2.isNext()) {
408  auto p = it2.get();
409  Point<3, double> xp = Particles.getPos(p);
410  Point<3, double> xP = Particles.getProp<8>(p);
411  Particles.getProp<0>(p) =0;
412  if (xP[0]==1.0) {
413 // Surface.add();
414 // Surface.last().get<0>() = p.getKey();
415  Particles.getProp<0>(p) = 0;
416  std::vector<double> SVel;
417  SVel=openfpm::math::sumY<K>(xP[0],xP[1],xP[2],Vr,V1,V2);
418  double SP=openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Vr);
419  Particles.getProp<2>(p)[0] = SVel[0];
420  Particles.getProp<2>(p)[1] = SVel[1];
421  Particles.getProp<2>(p)[2] = SVel[2];
422  Particles.getProp<9>(p)[0] = SVel[0];
423  Particles.getProp<9>(p)[1] = SVel[1];
424  Particles.getProp<9>(p)[2] = SVel[2];
425  Particles.getProp<5>(p) = SP;
426  Particles.setSubset(p,1);
427 
428  }
429  else {
430 // bulk.add();
431 // bulk.last().get<0>() = p.getKey();
432  Particles.setSubset(p,0);
433  Particles.getProp<0>(p) = 0;
434  Particles.getProp<1>(p)[0] = 0;
435  Particles.getProp<1>(p)[1] = 0;
436  Particles.getProp<1>(p)[2] = 0;
437  }
438  ++it2;
439  }
440 
441  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);
442  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);
443 
444  auto & bulk = Particles_bulk.getIds();
445  auto & Surface = Particles_surface.getIds();
446 
447  for (int j = 0; j < bulk.size(); j++) {
448  auto p = bulk.get<0>(j);
449  Point<3, double> xp = Particles.getPos(p);
450  Point<3, double> xP = Particles.getProp<8>(p);
451 
452  std::unordered_map<const lm,double,key_hash,key_equal> Ur;
453  std::unordered_map<const lm,double,key_hash,key_equal> U2;
454  std::unordered_map<const lm,double,key_hash,key_equal> U1;
455  std::unordered_map<const lm,double,key_hash,key_equal> Plm;
456 
457  for (int l = 0; l <= K; l++) {
458  for (int m = -l; m <= l; m++) {
459  auto Er= Vr.find(std::make_tuple(l,m));
460  auto E1= V1.find(std::make_tuple(l,m));
461  auto E2= V2.find(std::make_tuple(l,m));
462  std::vector<double> Sol=openfpm::math::sph_anasol_u(nu,l,m,Er->second,E1->second,E2->second,xP[0]);
463  Ur[std::make_tuple(l,m)]=Sol[0];
464  U1[std::make_tuple(l,m)]=Sol[1];
465  U2[std::make_tuple(l,m)]=Sol[2];
466  Plm[std::make_tuple(l,m)]=Sol[3];
467  }
468 
469  }
470 
471  if(fabs(xP[0])>=1e-5 && xP[1]>1e-5 && (M_PI-xP[1])>=1e-5)
472  {
473  std::vector<double> SVel = openfpm::math::sumY<K>(xP[0], xP[1], xP[2], Ur, U1, U2);
474  Particles.getProp<9>(p)[0] = SVel[0];
475  Particles.getProp<9>(p)[1] = SVel[1];
476  Particles.getProp<9>(p)[2] = SVel[2];
477  Particles.getProp<5>(p) = openfpm::math::sumY_Scalar<K>(xP[0], xP[1], xP[2], Plm);
478  }
479  }
480 
481  auto P = getV<0>(Particles);
482  auto V = getV<1>(Particles);
483  auto V_B = getV<2>(Particles);
484  V.setVarId(0);
485  auto DIV = getV<3>(Particles);
486  auto V_t = getV<4>(Particles);
487  auto P_anal = getV<5>(Particles);
488  auto temp=getV<6>(Particles);
489  auto RHS = getV<7>(Particles);
490  auto P_bulk = getV<0>(Particles_bulk);
491  auto RHS_bulk = getV<7>(Particles_bulk);
492  auto V_anal = getV<9>(Particles);
493 
494  V_t=V;
495  P=0;
496  P_bulk=0;
497  eq_id vx,vy,vz;
498 
499  vx.setId(0);
500  vy.setId(1);
501  vz.setId(2);
502 
503  double sampling=3.1;
504  double sampling2=1.9;
505  double rCut2=3.9*spacing;
506 
507  Particles.ghost_get_subset();
508  auto verletList = Particles.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
509  auto verletList2 = Particles.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut2);
510  auto verletList_bulk = Particles_bulk.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
511 
512  Derivative_x<decltype(verletList)> Dx(Particles, verletList, 2, rCut, support_options::RADIUS);
513  Derivative_y<decltype(verletList)> Dy(Particles, verletList, 2, rCut, support_options::RADIUS);
514  Derivative_z<decltype(verletList)> Dz(Particles, verletList, 2, rCut, support_options::RADIUS);
515 
516  Derivative_x<decltype(verletList_bulk)> B_Dx(Particles_bulk, verletList_bulk, 2, rCut, support_options::RADIUS);
517  Derivative_y<decltype(verletList_bulk)> B_Dy(Particles_bulk, verletList_bulk, 2, rCut, support_options::RADIUS);
518  Derivative_z<decltype(verletList_bulk)> B_Dz(Particles_bulk, verletList_bulk, 2, rCut, support_options::RADIUS);
519 
520  Derivative_xx<decltype(verletList2)> Dxx(Particles, verletList2, 2, rCut2, support_options::RADIUS);
521  Derivative_yy<decltype(verletList2)> Dyy(Particles, verletList2, 2, rCut2, support_options::RADIUS);
522  Derivative_zz<decltype(verletList2)> Dzz(Particles, verletList2, 2, rCut2, support_options::RADIUS);
523 
524  //std::cout << "DCPSE KERNELS DONE" << std::endl;
525  petsc_solver<double> solverPetsc;
526  solverPetsc.setPreconditioner(PCNONE);
527  timer tt;
528  double sum=0,sum1=0;
529  V_t=V;
530  double V_err_eps = 1e-5;
531 
532  double V_err = 1, V_err_old;
533  int n = 0;
534  int nmax = 30;
535  int ctr = 0, errctr, Vreset = 0;
536  V_err = 1;
537  n = 0;
538  double solvetime=0;
539  while (V_err >= V_err_eps && n <= nmax) {
540  //Particles.write_frame("StokesSphere",n);
541  Particles.ghost_get<0>(SKIP_LABELLING);
542  RHS_bulk[0] = B_Dx(P);
543  RHS_bulk[1] = B_Dy(P);
544  RHS_bulk[2] = B_Dz(P);
545  DCPSE_scheme<equations3d3, decltype(Particles)> Solver(Particles);
546  auto Stokes1 = nu * (Dxx(V[0])+Dyy(V[0])+Dzz(V[0]));
547  auto Stokes2 = nu * (Dxx(V[1])+Dyy(V[1])+Dzz(V[1]));
548  auto Stokes3 = nu * (Dxx(V[2])+Dyy(V[2])+Dzz(V[2]));
549  Solver.impose(Stokes1, bulk, RHS[0], vx);
550  Solver.impose(Stokes2, bulk, RHS[1], vy);
551  Solver.impose(Stokes3, bulk, RHS[2], vz);
552  Solver.impose(V[0], Surface, V_B[0], vx);
553  Solver.impose(V[1], Surface, V_B[1], vy);
554  Solver.impose(V[2], Surface, V_B[2], vz);
555  tt.start();
556  Solver.solve_with_solver(solverPetsc, V[0], V[1], V[2]);
557  tt.stop();
558  solvetime+=tt.getwct();
559  //Solver.solve(V[0],V[1],V[2]);
560  //std::cout << "Stokes Solved" << std::endl;
561  Particles.ghost_get<1>();
562  DIV = -(Dx(V[0])+Dy(V[1])+Dz(V[2]));
563  P_bulk = P + DIV;
564  sum = 0;
565  sum1 = 0;
566  for (int j = 0; j < bulk.size(); j++) {
567  auto p = bulk.get<0>(j);
568  sum += (Particles.getProp<4>(p)[0] - Particles.getProp<1>(p)[0]) *
569  (Particles.getProp<4>(p)[0] - Particles.getProp<1>(p)[0]) +
570  (Particles.getProp<4>(p)[1] - Particles.getProp<1>(p)[1]) *
571  (Particles.getProp<4>(p)[1] - Particles.getProp<1>(p)[1]) +
572  (Particles.getProp<4>(p)[2] - Particles.getProp<1>(p)[2]) *
573  (Particles.getProp<4>(p)[2] - Particles.getProp<1>(p)[2]);
574  sum1 += Particles.getProp<1>(p)[0] * Particles.getProp<1>(p)[0] +
575  Particles.getProp<1>(p)[1] * Particles.getProp<1>(p)[1] +
576  Particles.getProp<1>(p)[2] * Particles.getProp<1>(p)[2];
577  }
578  sum = sqrt(sum);
579  sum1 = sqrt(sum1);
580  v_cl.sum(sum);
581  v_cl.sum(sum1);
582  v_cl.execute();
583  V_t = V;
584  Particles.ghost_get<1>(SKIP_LABELLING);
585  V_err_old = V_err;
586  V_err = sum / sum1;
587  if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-14) {
588  errctr++;
589  } else {
590  errctr = 0;
591  }
592  if (n > 3) {
593  if (errctr > 1) {
594  std::cout << "CONVERGENCE LOOP BROKEN DUE TO INCREASE/VERY SLOW DECREASE IN ERROR" << std::endl;
595  Vreset = 1;
596  break;
597  } else {
598  Vreset = 0;
599  }
600  }
601  n++;
602 
603  }
604  //std::cout << "Total Solver time (wct):"<<solvetime<< std::endl;
605  V_t=0;
606 
607  double worst=0;
608  double L2=0;
609  for (int j = 0; j < bulk.size(); j++) {
610  auto p = bulk.get<0>(j);
611  Point<3,double> xP=Particles.getProp<8>(p);
612  if(xP[0]>=1e-5 && xP[1]>1e-5 && (M_PI-xP[1])>=1e-5)
613  {
614  double dx=Particles.getProp<9>(p)[0] - Particles.getProp<1>(p)[0];
615  double dy=Particles.getProp<9>(p)[1] - Particles.getProp<1>(p)[1];
616  double dz=Particles.getProp<9>(p)[2] - Particles.getProp<1>(p)[2];
617  Particles.getProp<4>(p)[0]=fabs(dx);
618  Particles.getProp<4>(p)[1]=fabs(dy);
619  Particles.getProp<4>(p)[2]=fabs(dz);
620  L2 += dx*dx+dy*dy+dz*dz;
621  if (std::max({fabs(dx),fabs(dy),fabs(dz)}) > worst) {
622  worst = std::max({fabs(dx),fabs(dy),fabs(dz)});
623  }
624 
625  }
626  }
627 
628  v_cl.sum(worst);
629  v_cl.sum(L2);
630  v_cl.execute();
631 /*
632  if (v_cl.rank() == 0) {
633  std::cout<<"Gd,Surf,Bulk Size: "<<grd_sz<<","<<Surface.size()<<","<<bulk.size()<<std::endl;
634  std::cout << "L2_Final: " <<sqrt(L2)<<","<<sqrt(L2/(bulk.size()+Surface.size()))
635  << std::endl;
636  std::cout << "L_inf_Final: " << worst
637  << std::endl;
638  }
639  std::cout << "L_inf_Final_test: " << worst;
640  Particles.write("StokesSphere");*/
641  BOOST_REQUIRE(worst<1e-3);
642 
643  }
644 
645 
646  BOOST_AUTO_TEST_CASE(Sph_harm_ig) {
647  BOOST_REQUIRE(openfpm::math::Y(2,1,0.5,0)+0.459674<0.00001);
648  //These would be a requirement once Boost releases their fix
649  //
650  //BOOST_REQUIRE(boost::math::legendre_p(0,-1,1)=?);
651  double nu=1.0;
652  size_t grd_sz=13;
653  const size_t sz[3] = {grd_sz,grd_sz,grd_sz};
654  Box<3, double> box({-1.0, -1.0,-1.0}, {1.0,1.0,1.0});
655  size_t bc[3] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
656  double spacing = 2.0 / (sz[0] - 1);
657  double rCut = 3.9*spacing;
658  double R=1.0;
659  Ghost<3, double> ghost(rCut);
660  // P V v_B RHS V_t P_anal RHS2 Polar cord
661  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);
662 
663 
664  auto &v_cl = create_vcluster();
665 
666 // openfpm::vector<aggregate<int>> bulk;
667 // openfpm::vector<aggregate<int>> Surface;
668 
669  auto it = Particles.getGridIterator(sz);
670  while (it.isNext()) {
671  auto key = it.get();
672  double x = -1.0+key.get(0) * it.getSpacing(0);
673  double y = -1.0+key.get(1) * it.getSpacing(1);
674  double z = -1.0+key.get(2) * it.getSpacing(2);
675  double r=sqrt(x*x+y*y+z*z);
676  if (r<R-spacing/2.0) {
677  Particles.add();
678  Particles.getLastPos()[0] = x;
679  Particles.getLastPos()[1] = y;
680  Particles.getLastPos()[2] = z;
681  Particles.getLastProp<8>()[0] = r;
682  if (r==0){
683  Particles.getLastProp<8>()[1] = 0.0;
684  }
685  else{
686  Particles.getLastProp<8>()[1] = std::atan2(sqrt(x*x+y*y),z);
687  }
688  Particles.getLastProp<8>()[2] = std::atan2(y,x);
689  }
690  ++it;
691  }
692 
693  int n_sp=int(grd_sz)*int(grd_sz)*3;
694 
695  double Golden_angle=M_PI * (3.0 - sqrt(5.0));
696 
697  for(int i=1;i<=n_sp;i++)
698  {
699  double y = 1.0 - (i /double(n_sp - 1.0)) * 2.0;
700  double radius = sqrt(1 - y * y);
701  double Golden_theta = Golden_angle * i;
702  double x = cos(Golden_theta) * radius;
703  double z = sin(Golden_theta) * radius;
704 
705  if (acos(z)==0 || acos(z)==M_PI){
706  std::cout<<"Theta 0/Pi "<<std::endl;
707  continue;
708  }
709 
710  Particles.add();
711  Particles.getLastPos()[0] = x;
712  Particles.getLastPos()[1] = y;
713  Particles.getLastPos()[2] = z;
714  Particles.getLastProp<8>()[0] = 1.0 ;
715  Particles.getLastProp<8>()[1] = std::atan2(sqrt(x*x+y*y),z);
716  Particles.getLastProp<8>()[2] = std::atan2(y,x);
717  }
718  Particles.map();
719  Particles.ghost_get<0>();
720 
721 
722  std::unordered_map<const lm,double,key_hash,key_equal> Vr;
723  std::unordered_map<const lm,double,key_hash,key_equal> V1;
724  std::unordered_map<const lm,double,key_hash,key_equal> V2;
725  //Setting max mode l_max
726  constexpr int K = 2;
727  //Setting amplitudes to 0
728  for(int l=0;l<=K;l++){
729  for(int m=-l;m<=l;m++){
730  Vr[std::make_tuple(l,m)]=0.0;
731  V1[std::make_tuple(l,m)]=0.0;
732  V2[std::make_tuple(l,m)]=0.0;
733  }
734 
735 
736  }
737  //Setting some amplitude for boundary velocity
738  V1[std::make_tuple(1,0)]=1.0;
739 
740  auto it2 = Particles.getDomainIterator();
741  while (it2.isNext()) {
742  auto p = it2.get();
743  Point<3, double> xp = Particles.getPos(p);
744  Point<3, double> xP = Particles.getProp<8>(p);
745  Particles.getProp<0>(p) =0;
746  if (xP[0]==1.0) {
747 // Surface.add();
748 // Surface.last().get<0>() = p.getKey();
749  Particles.getProp<0>(p) = 0;
750  std::vector<double> SVel;
751  SVel=openfpm::math::sumY<K>(xP[0],xP[1],xP[2],Vr,V1,V2);
752  double SP=openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Vr);
753  Particles.getProp<2>(p)[0] = SVel[0];
754  Particles.getProp<2>(p)[1] = SVel[1];
755  Particles.getProp<2>(p)[2] = SVel[2];
756  Particles.getProp<9>(p)[0] = SVel[0];
757  Particles.getProp<9>(p)[1] = SVel[1];
758  Particles.getProp<9>(p)[2] = SVel[2];
759  Particles.getProp<5>(p) = SP;
760  Particles.setSubset(p,1);
761 
762  }
763  else {
764 // bulk.add();
765 // bulk.last().get<0>() = p.getKey();
766  Particles.setSubset(p,0);
767  Particles.getProp<0>(p) = 0;
768  Particles.getProp<1>(p)[0] = 0;
769  Particles.getProp<1>(p)[1] = 0;
770  Particles.getProp<1>(p)[2] = 0;
771  }
772  ++it2;
773  }
774 
775  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);
776  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);
777 
778  auto & bulk = Particles_bulk.getIds();
779  auto & Surface = Particles_surface.getIds();
780 
781  for (int j = 0; j < bulk.size(); j++) {
782  auto p = bulk.get<0>(j);
783  Point<3, double> xp = Particles.getPos(p);
784  Point<3, double> xP = Particles.getProp<8>(p);
785 
786  std::unordered_map<const lm,double,key_hash,key_equal> Ur;
787  std::unordered_map<const lm,double,key_hash,key_equal> U2;
788  std::unordered_map<const lm,double,key_hash,key_equal> U1;
789  std::unordered_map<const lm,double,key_hash,key_equal> Plm;
790 
791  for (int l = 0; l <= K; l++) {
792  for (int m = -l; m <= l; m++) {
793  auto Er= Vr.find(std::make_tuple(l,m));
794  auto E1= V1.find(std::make_tuple(l,m));
795  auto E2= V2.find(std::make_tuple(l,m));
796  std::vector<double> Sol=openfpm::math::sph_anasol_u(nu,l,m,Er->second,E1->second,E2->second,xP[0]);
797  Ur[std::make_tuple(l,m)]=Sol[0];
798  U1[std::make_tuple(l,m)]=Sol[1];
799  U2[std::make_tuple(l,m)]=Sol[2];
800  Plm[std::make_tuple(l,m)]=Sol[3];
801  }
802 
803  }
804 
805  if(fabs(xP[0])>=1e-5 && xP[1]>1e-5 && (M_PI-xP[1])>=1e-5)
806  {
807  std::vector<double> SVel = openfpm::math::sumY<K>(xP[0], xP[1], xP[2], Ur, U1, U2);
808  Particles.getProp<9>(p)[0] = SVel[0];
809  Particles.getProp<9>(p)[1] = SVel[1];
810  Particles.getProp<9>(p)[2] = SVel[2];
811  Particles.getProp<5>(p) = openfpm::math::sumY_Scalar<K>(xP[0], xP[1], xP[2], Plm);
812  }
813  }
814 
815  auto P = getV<0>(Particles);
816  auto V = getV<1>(Particles);
817  auto V_B = getV<2>(Particles);
818  V.setVarId(0);
819  auto DIV = getV<3>(Particles);
820  auto V_t = getV<4>(Particles);
821  auto P_anal = getV<5>(Particles);
822  auto temp=getV<6>(Particles);
823  auto RHS = getV<7>(Particles);
824  auto P_bulk = getV<0>(Particles_bulk);
825  auto RHS_bulk = getV<7>(Particles_bulk);
826  auto V_anal = getV<9>(Particles);
827 
828  V_t=V;
829  P=0;
830  P_bulk=0;
831  eq_id vx,vy,vz;
832 
833  vx.setId(0);
834  vy.setId(1);
835  vz.setId(2);
836 
837  double sampling=3.1;
838  double sampling2=1.9;
839  double rCut2=3.9*spacing;
840 
841  Particles.ghost_get_subset();
842  auto verletList = Particles.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
843  auto verletList2 = Particles.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut2);
844  auto verletList_bulk = Particles_bulk.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut2);
845 
846  Derivative_x<decltype(verletList)> Dx(Particles, verletList, 2, rCut, support_options::RADIUS);
847  Derivative_y<decltype(verletList)> Dy(Particles, verletList, 2, rCut, support_options::RADIUS);
848  Derivative_z<decltype(verletList)> Dz(Particles, verletList, 2, rCut, support_options::RADIUS);
849 
850  Derivative_x<decltype(verletList_bulk)> B_Dx(Particles_bulk, verletList_bulk, 2, rCut, support_options::RADIUS);
851  Derivative_y<decltype(verletList_bulk)> B_Dy(Particles_bulk, verletList_bulk, 2, rCut, support_options::RADIUS);
852  Derivative_z<decltype(verletList_bulk)> B_Dz(Particles_bulk, verletList_bulk, 2, rCut, support_options::RADIUS);
853 
854  Derivative_xx<decltype(verletList2)> Dxx(Particles, verletList2, 2, rCut2, support_options::RADIUS);
855  Derivative_yy<decltype(verletList2)> Dyy(Particles, verletList2, 2, rCut2, support_options::RADIUS);
856  Derivative_zz<decltype(verletList2)> Dzz(Particles, verletList2, 2, rCut2, support_options::RADIUS);
857 
858  //std::cout << "DCPSE KERNELS DONE" << std::endl;
859  petsc_solver<double> solverPetsc;
860  solverPetsc.setPreconditioner(PCNONE);
861  timer tt;
862  double sum=0,sum1=0;
863  V_t=V;
864  double V_err_eps = 1e-5;
865 
866  double V_err = 1, V_err_old;
867  int n = 0;
868  int nmax = 30;
869  int ctr = 0, errctr, Vreset = 0;
870  V_err = 1;
871  n = 0;
872  double solvetime=0;
873  while (V_err >= V_err_eps && n <= nmax) {
874  //Particles.write_frame("StokesSphere",n);
875  Particles.ghost_get<0>(SKIP_LABELLING);
876  RHS_bulk[0] = B_Dx(P);
877  RHS_bulk[1] = B_Dy(P);
878  RHS_bulk[2] = B_Dz(P);
879  DCPSE_scheme<equations3d3, decltype(Particles)> Solver(Particles);
880  auto Stokes1 = nu * (Dxx(V[0])+Dyy(V[0])+Dzz(V[0]));
881  auto Stokes2 = nu * (Dxx(V[1])+Dyy(V[1])+Dzz(V[1]));
882  auto Stokes3 = nu * (Dxx(V[2])+Dyy(V[2])+Dzz(V[2]));
883  Solver.impose(Stokes1, bulk, RHS[0], vx);
884  Solver.impose(Stokes2, bulk, RHS[1], vy);
885  Solver.impose(Stokes3, bulk, RHS[2], vz);
886  Solver.impose(V[0], Surface, V_B[0], vx);
887  Solver.impose(V[1], Surface, V_B[1], vy);
888  Solver.impose(V[2], Surface, V_B[2], vz);
889  Solver.impose_x_ig(bulk, V[0], vx);
890  Solver.impose_x_ig(bulk, V[1], vy);
891  Solver.impose_x_ig(bulk, V[2], vz);
892  Solver.impose_x_ig(Surface, V[0], vx);
893  Solver.impose_x_ig(Surface, V[1], vy);
894  Solver.impose_x_ig(Surface, V[2], vz);
895  tt.start();
896  Solver.solve_with_solver_ig(solverPetsc, V[0], V[1], V[2]);
897  tt.stop();
898  solvetime+=tt.getwct();
899  //Solver.solve(V[0],V[1],V[2]);
900  //std::cout << "Stokes Solved" << std::endl;
901  Particles.ghost_get<1>();
902  DIV = -(Dx(V[0])+Dy(V[1])+Dz(V[2]));
903  P_bulk = P + DIV;
904  sum = 0;
905  sum1 = 0;
906  for (int j = 0; j < bulk.size(); j++) {
907  auto p = bulk.get<0>(j);
908  sum += (Particles.getProp<4>(p)[0] - Particles.getProp<1>(p)[0]) *
909  (Particles.getProp<4>(p)[0] - Particles.getProp<1>(p)[0]) +
910  (Particles.getProp<4>(p)[1] - Particles.getProp<1>(p)[1]) *
911  (Particles.getProp<4>(p)[1] - Particles.getProp<1>(p)[1]) +
912  (Particles.getProp<4>(p)[2] - Particles.getProp<1>(p)[2]) *
913  (Particles.getProp<4>(p)[2] - Particles.getProp<1>(p)[2]);
914  sum1 += Particles.getProp<1>(p)[0] * Particles.getProp<1>(p)[0] +
915  Particles.getProp<1>(p)[1] * Particles.getProp<1>(p)[1] +
916  Particles.getProp<1>(p)[2] * Particles.getProp<1>(p)[2];
917  }
918  sum = sqrt(sum);
919  sum1 = sqrt(sum1);
920  v_cl.sum(sum);
921  v_cl.sum(sum1);
922  v_cl.execute();
923  V_t = V;
924  Particles.ghost_get<1>(SKIP_LABELLING);
925  V_err_old = V_err;
926  V_err = sum / sum1;
927  if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-14) {
928  errctr++;
929  } else {
930  errctr = 0;
931  }
932  if (n > 3) {
933  if (errctr > 1) {
934  std::cout << "CONVERGENCE LOOP BROKEN DUE TO INCREASE/VERY SLOW DECREASE IN ERROR" << std::endl;
935  Vreset = 1;
936  break;
937  } else {
938  Vreset = 0;
939  }
940  }
941  n++;
942 
943  }
944  //std::cout << "Total Solver time (wct):"<<solvetime<< std::endl;
945 
946  V_t=0;
947 
948  double worst=0;
949  double L2=0;
950  for (int j = 0; j < bulk.size(); j++) {
951  auto p = bulk.get<0>(j);
952  Point<3,double> xP=Particles.getProp<8>(p);
953  if(xP[0]>=1e-5 && xP[1]>1e-5 && (M_PI-xP[1])>=1e-5)
954  {
955  double dx=Particles.getProp<9>(p)[0] - Particles.getProp<1>(p)[0];
956  double dy=Particles.getProp<9>(p)[1] - Particles.getProp<1>(p)[1];
957  double dz=Particles.getProp<9>(p)[2] - Particles.getProp<1>(p)[2];
958  Particles.getProp<4>(p)[0]=fabs(dx);
959  Particles.getProp<4>(p)[1]=fabs(dy);
960  Particles.getProp<4>(p)[2]=fabs(dz);
961  L2 += dx*dx+dy*dy+dz*dz;
962  if (std::max({fabs(dx),fabs(dy),fabs(dz)}) > worst) {
963  worst = std::max({fabs(dx),fabs(dy),fabs(dz)});
964  }
965 
966  }
967  }
968 
969  v_cl.sum(worst);
970  v_cl.sum(L2);
971  v_cl.execute();
972 
973  /* if (v_cl.rank() == 0) {
974  std::cout<<"Gd,Surf,Bulk Size: "<<grd_sz<<","<<Surface.size()<<","<<bulk.size()<<std::endl;
975  std::cout << "L2_Final: " <<sqrt(L2)<<","<<sqrt(L2/(bulk.size()+Surface.size()))
976  << std::endl;
977  std::cout << "L_inf_Final: " << worst
978  << std::endl;
979  }
980  std::cout << "L_inf_Final_test: " << worst;*/
981  //Particles.write("StokesSphere");
982  BOOST_REQUIRE(worst<1e-3);
983 
984  }
985 BOOST_AUTO_TEST_SUITE_END()
986 
987 #endif
988 #endif
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_test
Test structure used for several test.
Definition: Point_test.hpp:106
Point
This class implement the point shape in an N-dimensional space.
Definition: Point.hpp:28
Point::get
__device__ __host__ const T & get(unsigned int i) const
Get coordinate.
Definition: Point.hpp:172
eq_id
Definition: eq_solve_common.hpp:15
openfpm::vector
Implementation of 1-D std::vector like structure.
Definition: map_vector.hpp: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
timer
Class for cpu time benchmarking.
Definition: timer.hpp:28
timer::stop
void stop()
Stop the timer.
Definition: timer.hpp:119
timer::start
void start()
Start the timer.
Definition: timer.hpp:90
timer::getwct
double getwct()
Return the elapsed real time.
Definition: timer.hpp:130
vector_dist_subset
Definition: vector_dist_subset.hpp:84
vector_dist_ws
Definition: vector_dist_subset.hpp:17
vector_dist
Distributed vector.
Definition: vector_dist.hpp:176
cub::int
KeyT const ValueT ValueT OffsetIteratorT OffsetIteratorT int
[in] The number of segments that comprise the sorting data
Definition: dispatch_radix_sort.cuh:336
prop_id
Definition: eq_solve_common.hpp:12
sum
It model an expression expr1 + ... exprn.
Definition: sum.hpp:93