OpenFPM  5.2.0
Project that contain the implementation of distributed structures
main.cpp
1 //
2 // Created by Abhinav Singh on 15.06.20.
3 //
4 
40 // Include Vector Expression,Vector Expressions for Subset,DCPSE , and Solver header files
41 #include "config.h"
42 #include <iostream>
43 #include "DCPSE/DCPSE_op/DCPSE_op.hpp"
44 #include "DCPSE/DCPSE_op/DCPSE_Solver.hpp"
45 #include "Operators/Vector/vector_dist_operators.hpp"
46 #include "Vector/vector_dist_subset.hpp"
47 #include "util/SphericalHarmonics.hpp"
49 
50 int main(int argc, char* argv[])
51 {
52 
53  {
70  openfpm_init(&argc,&argv);
71  timer tt2;
72  tt2.start();
73  size_t grd_sz = 18;
74  double V_err_eps = 5e-4;
75 
76  double nu=1.0;
77  const size_t sz[3] = {grd_sz,grd_sz,grd_sz};
78  Box<3, double> box({-1.0, -1.0,-1.0}, {1.0,1.0,1.0});
79  size_t bc[3] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
80  double spacing = 2.0 / (sz[0] - 1);
81  double rCut = 3.9*spacing;
82  double R=1.0;
83  Ghost<3, double> ghost(rCut);
84  // P V v_B RHS V_t P_anal RHS2 Polar cord
85  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);
86  Particles.setPropNames({"Pressure","Velocity","Velocity_Boundary","Divergence","V_T","P_analytical","TEMP" ,"RHS","PolarCoord","V_analytical"});
88 
89 
101  auto &v_cl = create_vcluster();
102 
103  auto it = Particles.getGridIterator(sz);
104  while (it.isNext()) {
105  auto key = it.get();
106  double x = -1.0+key.get(0) * it.getSpacing(0);
107  double y = -1.0+key.get(1) * it.getSpacing(1);
108  double z = -1.0+key.get(2) * it.getSpacing(2);
109  double r=sqrt(x*x+y*y+z*z);
110  if (r<R-spacing/2.0) {
111  Particles.add();
112  Particles.getLastPos()[0] = x;
113  Particles.getLastPos()[1] = y;
114  Particles.getLastPos()[2] = z;
115  Particles.getLastProp<8>()[0] = r;
116  if (r==0){
117  Particles.getLastProp<8>()[1] = 0.0;
118  }
119  else{
120  Particles.getLastProp<8>()[1] = std::atan2(sqrt(x*x+y*y),z);
121  }
122  Particles.getLastProp<8>()[2] = std::atan2(y,x);
123  }
124  ++it;
125  }
126 
127  int n_sp=int(grd_sz)*int(grd_sz)*3;
128 
129  double Golden_angle=M_PI * (3.0 - sqrt(5.0));
130 
131  for(int i=1;i<=n_sp;i++)
132  {
133  double y = 1.0 - (i /double(n_sp - 1.0)) * 2.0;
134  double radius = sqrt(1 - y * y);
135  double Golden_theta = Golden_angle * i;
136  double x = cos(Golden_theta) * radius;
137  double z = sin(Golden_theta) * radius;
138 
139  if (acos(z)==0 || acos(z)==M_PI){
140  std::cout<<"Theta 0/Pi "<<std::endl;
141  continue;
142  }
143 
144  Particles.add();
145  Particles.getLastPos()[0] = x;
146  Particles.getLastPos()[1] = y;
147  Particles.getLastPos()[2] = z;
148  Particles.getLastProp<8>()[0] = 1.0 ;
149  Particles.getLastProp<8>()[1] = std::atan2(sqrt(x*x+y*y),z);
150  Particles.getLastProp<8>()[2] = std::atan2(y,x);
151  }
152  Particles.map();
153  Particles.ghost_get<0>();
155 
167 
168  std::unordered_map<const lm,double,key_hash,key_equal> Vr;
169  std::unordered_map<const lm,double,key_hash,key_equal> V1;
170  std::unordered_map<const lm,double,key_hash,key_equal> V2;
171  //Setting max mode l_max
172  constexpr int K = 2;
173  //Setting amplitudes to 0
174  for(int l=0;l<=K;l++){
175  for(int m=-l;m<=l;m++){
176  Vr[std::make_tuple(l,m)]=0.0;
177  V1[std::make_tuple(l,m)]=0.0;
178  V2[std::make_tuple(l,m)]=0.0;
179  }
180 
181 
182  }
183  //Setting some amplitude for boundary velocity
184  V1[std::make_tuple(2,0)]=1.0;
185 
186  auto it2 = Particles.getDomainIterator();
187  while (it2.isNext()) {
188  auto p = it2.get();
189  Point<3, double> xp = Particles.getPos(p);
190  Point<3, double> xP = Particles.getProp<8>(p);
191  Particles.getProp<0>(p) =0;
192  if (xP[0]==1.0) {
193  Particles.getProp<0>(p) = 0;
194  std::vector<double> SVel;
195  SVel=openfpm::math::sumY<K>(xP[0],xP[1],xP[2],Vr,V1,V2);
196  double SP=openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Vr);
197  Particles.getProp<2>(p)[0] = SVel[0];
198  Particles.getProp<2>(p)[1] = SVel[1];
199  Particles.getProp<2>(p)[2] = SVel[2];
200  Particles.getProp<9>(p)[0] = SVel[0];
201  Particles.getProp<9>(p)[1] = SVel[1];
202  Particles.getProp<9>(p)[2] = SVel[2];
203  Particles.getProp<5>(p) = SP;
204  Particles.setSubset(p,1);
205 
206  }
207  else {
208  Particles.setSubset(p,0);
209  Particles.getProp<0>(p) = 0;
210  Particles.getProp<1>(p)[0] = 0;
211  Particles.getProp<1>(p)[1] = 0;
212  Particles.getProp<1>(p)[2] = 0;
213  }
214  ++it2;
215  }
217 
228 
229  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);
230  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);
231 
232  auto & bulk = Particles_bulk.getIds();
233  auto & Surface = Particles_surface.getIds();
234 
235  for (int j = 0; j < bulk.size(); j++) {
236  auto p = bulk.get<0>(j);
237  Point<3, double> xp = Particles.getPos(p);
238  Point<3, double> xP = Particles.getProp<8>(p);
239 
240  //Dictionary for basis coefficients
241  std::unordered_map<const lm,double,key_hash,key_equal> Ur, U1, U2, Plm;
242 
243  //Setting basis coefficients
244  for (int l = 0; l <= K; l++) {
245  for (int m = -l; m <= l; m++) {
246  auto Er= Vr.find(std::make_tuple(l,m));
247  auto E1= V1.find(std::make_tuple(l,m));
248  auto E2= V2.find(std::make_tuple(l,m));
249  std::vector<double> Sol=openfpm::math::sph_anasol_u(nu,l,m,Er->second,E1->second,E2->second,xP[0]);
250  Ur[std::make_tuple(l,m)]=Sol[0];
251  U1[std::make_tuple(l,m)]=Sol[1];
252  U2[std::make_tuple(l,m)]=Sol[2];
253  Plm[std::make_tuple(l,m)]=Sol[3];
254  }
255 
256  }
257 
258  //Computing Analyical Solution.
259  if(fabs(xP[0])>=1e-5 && xP[1]>1e-5 && (M_PI-xP[1])>=1e-5)
260  {
261  std::vector<double> SVel = openfpm::math::sumY<K>(xP[0], xP[1], xP[2], Ur, U1, U2);
262  Particles.getProp<9>(p)[0] = SVel[0];
263  Particles.getProp<9>(p)[1] = SVel[1];
264  Particles.getProp<9>(p)[2] = SVel[2];
265  Particles.getProp<5>(p) = openfpm::math::sumY_Scalar<K>(xP[0], xP[1], xP[2], Plm);
266  }
267  }
269 
280  auto P = getV<0>(Particles);
281  auto V = getV<1>(Particles);
282  auto V_B = getV<2>(Particles);
283  V.setVarId(0);
284  auto DIV = getV<3>(Particles);
285  auto V_t = getV<4>(Particles);
286  auto P_anal = getV<5>(Particles);
287  auto temp=getV<6>(Particles);
288  auto RHS = getV<7>(Particles);
289  auto P_bulk = getV<0>(Particles_bulk);
290  auto RHS_bulk = getV<7>(Particles_bulk);
291  auto V_anal = getV<9>(Particles);
292 
293  V_t=V;
294  P=0;
295  P_bulk=0;
296  double sampling=3.1;
297  double sampling2=1.9;
298  double rCut2=3.9*spacing;
299 
300  auto verletList = Particles.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
301  auto verletList2 = Particles.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut2);
302  auto verletList_bulk = Particles_bulk.template getVerlet<VL_NON_SYMMETRIC|VL_SKIP_REF_PART>(rCut);
303 
304  Derivative_x<decltype(verletList)> Dx(Particles, verletList, 2, rCut, support_options::RADIUS);
305  Derivative_y<decltype(verletList)> Dy(Particles, verletList, 2, rCut, support_options::RADIUS);
306  Derivative_z<decltype(verletList)> Dz(Particles, verletList, 2, rCut, support_options::RADIUS);
307  Derivative_x<decltype(verletList_bulk)> B_Dx(Particles, Particles_bulk, verletList_bulk, 2, rCut, support_options::RADIUS);
308  Derivative_y<decltype(verletList_bulk)> B_Dy(Particles, Particles_bulk, verletList_bulk, 2, rCut, support_options::RADIUS);
309  Derivative_z<decltype(verletList_bulk)> B_Dz(Particles, Particles_bulk, verletList_bulk, 2, rCut, support_options::RADIUS);
310  Derivative_xx<decltype(verletList2)> Dxx(Particles, verletList2, 2, rCut2, support_options::RADIUS);
311  Derivative_yy<decltype(verletList2)> Dyy(Particles, verletList2, 2, rCut2, support_options::RADIUS);
312  Derivative_zz<decltype(verletList2)> Dzz(Particles, verletList2, 2, rCut2, support_options::RADIUS);
313 
315 
326 
327  eq_id vx,vy,vz;
328 
329  vx.setId(0);
330  vy.setId(1);
331  vz.setId(2);
332 
333  petsc_solver<double> solverPetsc;
334  solverPetsc.setPreconditioner(PCNONE);
335  timer tt;
336  double sum=0,sum1=0;
337  V_t=V;
338  double V_err = 1, V_err_old;
339  int n = 0, nmax = 30, ctr = 0, errctr, Vreset = 0;
340  V_err = 1;
341  n = 0;
342  tt.start();
343  while (V_err >= V_err_eps && n <= nmax) {
344  Particles.ghost_get<0>(SKIP_LABELLING);
345  RHS_bulk[0] = B_Dx(P);
346  RHS_bulk[1] = B_Dy(P);
347  RHS_bulk[2] = B_Dz(P);
348  DCPSE_scheme<equations3d3, decltype(Particles)> Solver(Particles);
349  auto Stokes1 = nu * (Dxx(V[0])+Dyy(V[0])+Dzz(V[0]));
350  auto Stokes2 = nu * (Dxx(V[1])+Dyy(V[1])+Dzz(V[1]));
351  auto Stokes3 = nu * (Dxx(V[2])+Dyy(V[2])+Dzz(V[2]));
352  Solver.impose(Stokes1, bulk, RHS[0], vx);
353  Solver.impose(Stokes2, bulk, RHS[1], vy);
354  Solver.impose(Stokes3, bulk, RHS[2], vz);
355  Solver.impose(V[0], Surface, V_B[0], vx);
356  Solver.impose(V[1], Surface, V_B[1], vy);
357  Solver.impose(V[2], Surface, V_B[2], vz);
358  Solver.solve_with_solver(solverPetsc, V[0], V[1], V[2]);
359  Particles.ghost_get<1>();
360  DIV = -(Dx(V[0])+Dy(V[1])+Dz(V[2]));
361  P_bulk = P + DIV;
362  sum = 0;
363  sum1 = 0;
364  for (int j = 0; j < bulk.size(); j++) {
365  auto p = bulk.get<0>(j);
366  sum += (Particles.getProp<4>(p)[0] - Particles.getProp<1>(p)[0]) *
367  (Particles.getProp<4>(p)[0] - Particles.getProp<1>(p)[0]) +
368  (Particles.getProp<4>(p)[1] - Particles.getProp<1>(p)[1]) *
369  (Particles.getProp<4>(p)[1] - Particles.getProp<1>(p)[1]) +
370  (Particles.getProp<4>(p)[2] - Particles.getProp<1>(p)[2]) *
371  (Particles.getProp<4>(p)[2] - Particles.getProp<1>(p)[2]);
372  sum1 += Particles.getProp<1>(p)[0] * Particles.getProp<1>(p)[0] +
373  Particles.getProp<1>(p)[1] * Particles.getProp<1>(p)[1] +
374  Particles.getProp<1>(p)[2] * Particles.getProp<1>(p)[2];
375  }
376  sum = sqrt(sum);
377  sum1 = sqrt(sum1);
378  v_cl.sum(sum);
379  v_cl.sum(sum1);
380  v_cl.execute();
381  V_t = V;
382  Particles.ghost_get<1>(SKIP_LABELLING);
383  V_err_old = V_err;
384  V_err = sum / sum1;
385  if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-14) {
386  errctr++;
387  } else {
388  errctr = 0;
389  }
390  if (n > 3) {
391  if (errctr > 1) {
392  std::cout << "CONVERGENCE LOOP BROKEN DUE TO INCREASE/VERY SLOW DECREASE IN Divergence" << std::endl;
393  Vreset = 1;
394  break;
395  } else {
396  Vreset = 0;
397  }
398  }
399  n++;
400 
401  }
402  //Writing the final Solution
403  //The solution can be visualized at https://link.springer.com/article/10.1140/epje/s10189-021-00121-x/figures/7
404  Particles.write("StokesSphere");
405  } //Ending Scope for Petsc.
406  //Finalizing the Library.
407  openfpm_finalize();
409 
418 }
Box
This class represent an N-dimensional box.
Definition: Box.hpp:60
Ghost
Definition: Ghost.hpp:40
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
eq_id
Definition: eq_solve_common.hpp:15
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::start
void start()
Start the timer.
Definition: timer.hpp:90
vector_dist_subset
Definition: vector_dist_subset.hpp:84
vector_dist_ws
Definition: vector_dist_subset.hpp:17
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
sum
It model an expression expr1 + ... exprn.
Definition: sum.hpp:93