OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
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
50int 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 Derivative_x Dx(Particles, 2, rCut,sampling, support_options::RADIUS),B_Dx(Particles_bulk, 2, rCut,sampling, support_options::RADIUS);
301 Derivative_y Dy(Particles, 2, rCut,sampling, support_options::RADIUS),B_Dy(Particles_bulk, 2, rCut,sampling, support_options::RADIUS);
302 Derivative_z Dz(Particles, 2, rCut,sampling, support_options::RADIUS),B_Dz(Particles_bulk, 2, rCut,sampling, support_options::RADIUS);
303 Derivative_xx Dxx(Particles, 2, rCut2,sampling2,support_options::RADIUS);
304 Derivative_yy Dyy(Particles, 2, rCut2,sampling2,support_options::RADIUS);
305 Derivative_zz Dzz(Particles, 2, rCut2,sampling2,support_options::RADIUS);
306
308
319
320 eq_id vx,vy,vz;
321
322 vx.setId(0);
323 vy.setId(1);
324 vz.setId(2);
325
326 petsc_solver<double> solverPetsc;
327 solverPetsc.setPreconditioner(PCNONE);
328 timer tt;
329 double sum=0,sum1=0;
330 V_t=V;
331 double V_err = 1, V_err_old;
332 int n = 0, nmax = 30, ctr = 0, errctr, Vreset = 0;
333 V_err = 1;
334 n = 0;
335 tt.start();
336 while (V_err >= V_err_eps && n <= nmax) {
337 Particles.ghost_get<0>(SKIP_LABELLING);
338 RHS_bulk[0] = B_Dx(P);
339 RHS_bulk[1] = B_Dy(P);
340 RHS_bulk[2] = B_Dz(P);
341 DCPSE_scheme<equations3d3, decltype(Particles)> Solver(Particles);
342 auto Stokes1 = nu * (Dxx(V[0])+Dyy(V[0])+Dzz(V[0]));
343 auto Stokes2 = nu * (Dxx(V[1])+Dyy(V[1])+Dzz(V[1]));
344 auto Stokes3 = nu * (Dxx(V[2])+Dyy(V[2])+Dzz(V[2]));
345 Solver.impose(Stokes1, bulk, RHS[0], vx);
346 Solver.impose(Stokes2, bulk, RHS[1], vy);
347 Solver.impose(Stokes3, bulk, RHS[2], vz);
348 Solver.impose(V[0], Surface, V_B[0], vx);
349 Solver.impose(V[1], Surface, V_B[1], vy);
350 Solver.impose(V[2], Surface, V_B[2], vz);
351 Solver.solve_with_solver(solverPetsc, V[0], V[1], V[2]);
352 Particles.ghost_get<1>();
353 DIV = -(Dx(V[0])+Dy(V[1])+Dz(V[2]));
354 P_bulk = P + DIV;
355 sum = 0;
356 sum1 = 0;
357 for (int j = 0; j < bulk.size(); j++) {
358 auto p = bulk.get<0>(j);
359 sum += (Particles.getProp<4>(p)[0] - Particles.getProp<1>(p)[0]) *
360 (Particles.getProp<4>(p)[0] - Particles.getProp<1>(p)[0]) +
361 (Particles.getProp<4>(p)[1] - Particles.getProp<1>(p)[1]) *
362 (Particles.getProp<4>(p)[1] - Particles.getProp<1>(p)[1]) +
363 (Particles.getProp<4>(p)[2] - Particles.getProp<1>(p)[2]) *
364 (Particles.getProp<4>(p)[2] - Particles.getProp<1>(p)[2]);
365 sum1 += Particles.getProp<1>(p)[0] * Particles.getProp<1>(p)[0] +
366 Particles.getProp<1>(p)[1] * Particles.getProp<1>(p)[1] +
367 Particles.getProp<1>(p)[2] * Particles.getProp<1>(p)[2];
368 }
369 sum = sqrt(sum);
370 sum1 = sqrt(sum1);
371 v_cl.sum(sum);
372 v_cl.sum(sum1);
373 v_cl.execute();
374 V_t = V;
375 Particles.ghost_get<1>(SKIP_LABELLING);
376 V_err_old = V_err;
377 V_err = sum / sum1;
378 if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-14) {
379 errctr++;
380 } else {
381 errctr = 0;
382 }
383 if (n > 3) {
384 if (errctr > 1) {
385 std::cout << "CONVERGENCE LOOP BROKEN DUE TO INCREASE/VERY SLOW DECREASE IN Divergence" << std::endl;
386 Vreset = 1;
387 break;
388 } else {
389 Vreset = 0;
390 }
391 }
392 n++;
393
394 }
395 //Writing the final Solution
396 //The solution can be visualized at https://link.springer.com/article/10.1140/epje/s10189-021-00121-x/figures/7
397 Particles.write("StokesSphere");
398 } //Ending Scope for Petsc.
399 //Finalizing the Library.
400 openfpm_finalize();
402
411}
Box
This class represent an N-dimensional box.
Definition Box.hpp:61
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
In case T does not match the PETSC precision compilation create a stub structure.
Definition petsc_solver.hpp:56
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