OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
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
32BOOST_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=13;
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 double solvetime=0;
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 tt.start();
544 Solver.solve_with_solver(solverPetsc, V[0], V[1], V[2]);
545 tt.stop();
546 solvetime+=tt.getwct();
547 //Solver.solve(V[0],V[1],V[2]);
548 //std::cout << "Stokes Solved" << std::endl;
549 Particles.ghost_get<1>();
550 DIV = -(Dx(V[0])+Dy(V[1])+Dz(V[2]));
551 P_bulk = P + DIV;
552 sum = 0;
553 sum1 = 0;
554 for (int j = 0; j < bulk.size(); j++) {
555 auto p = bulk.get<0>(j);
556 sum += (Particles.getProp<4>(p)[0] - Particles.getProp<1>(p)[0]) *
557 (Particles.getProp<4>(p)[0] - Particles.getProp<1>(p)[0]) +
558 (Particles.getProp<4>(p)[1] - Particles.getProp<1>(p)[1]) *
559 (Particles.getProp<4>(p)[1] - Particles.getProp<1>(p)[1]) +
560 (Particles.getProp<4>(p)[2] - Particles.getProp<1>(p)[2]) *
561 (Particles.getProp<4>(p)[2] - Particles.getProp<1>(p)[2]);
562 sum1 += Particles.getProp<1>(p)[0] * Particles.getProp<1>(p)[0] +
563 Particles.getProp<1>(p)[1] * Particles.getProp<1>(p)[1] +
564 Particles.getProp<1>(p)[2] * Particles.getProp<1>(p)[2];
565 }
566 sum = sqrt(sum);
567 sum1 = sqrt(sum1);
568 v_cl.sum(sum);
569 v_cl.sum(sum1);
570 v_cl.execute();
571 V_t = V;
572 Particles.ghost_get<1>(SKIP_LABELLING);
573 V_err_old = V_err;
574 V_err = sum / sum1;
575 if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-14) {
576 errctr++;
577 } else {
578 errctr = 0;
579 }
580 if (n > 3) {
581 if (errctr > 1) {
582 std::cout << "CONVERGENCE LOOP BROKEN DUE TO INCREASE/VERY SLOW DECREASE IN ERROR" << std::endl;
583 Vreset = 1;
584 break;
585 } else {
586 Vreset = 0;
587 }
588 }
589 n++;
590
591 }
592 //std::cout << "Total Solver time (wct):"<<solvetime<< std::endl;
593 V_t=0;
594
595 double worst=0;
596 double L2=0;
597 for (int j = 0; j < bulk.size(); j++) {
598 auto p = bulk.get<0>(j);
599 Point<3,double> xP=Particles.getProp<8>(p);
600 if(xP[0]>=1e-5 && xP[1]>1e-5 && (M_PI-xP[1])>=1e-5)
601 {
602 double dx=Particles.getProp<9>(p)[0] - Particles.getProp<1>(p)[0];
603 double dy=Particles.getProp<9>(p)[1] - Particles.getProp<1>(p)[1];
604 double dz=Particles.getProp<9>(p)[2] - Particles.getProp<1>(p)[2];
605 Particles.getProp<4>(p)[0]=fabs(dx);
606 Particles.getProp<4>(p)[1]=fabs(dy);
607 Particles.getProp<4>(p)[2]=fabs(dz);
608 L2 += dx*dx+dy*dy+dz*dz;
609 if (std::max({fabs(dx),fabs(dy),fabs(dz)}) > worst) {
610 worst = std::max({fabs(dx),fabs(dy),fabs(dz)});
611 }
612
613 }
614 }
615
616 v_cl.sum(worst);
617 v_cl.sum(L2);
618 v_cl.execute();
619/*
620 if (v_cl.rank() == 0) {
621 std::cout<<"Gd,Surf,Bulk Size: "<<grd_sz<<","<<Surface.size()<<","<<bulk.size()<<std::endl;
622 std::cout << "L2_Final: " <<sqrt(L2)<<","<<sqrt(L2/(bulk.size()+Surface.size()))
623 << std::endl;
624 std::cout << "L_inf_Final: " << worst
625 << std::endl;
626 }
627 std::cout << "L_inf_Final_test: " << worst;
628 Particles.write("StokesSphere");*/
629 BOOST_REQUIRE(worst<1e-3);
630
631 }
632
633
634 BOOST_AUTO_TEST_CASE(Sph_harm_ig) {
635 BOOST_REQUIRE(openfpm::math::Y(2,1,0.5,0)+0.459674<0.00001);
636 //These would be a requirement once Boost releases their fix
637 //
638 //BOOST_REQUIRE(boost::math::legendre_p(0,-1,1)=?);
639 double nu=1.0;
640 size_t grd_sz=13;
641 const size_t sz[3] = {grd_sz,grd_sz,grd_sz};
642 Box<3, double> box({-1.0, -1.0,-1.0}, {1.0,1.0,1.0});
643 size_t bc[3] = {NON_PERIODIC, NON_PERIODIC, NON_PERIODIC};
644 double spacing = 2.0 / (sz[0] - 1);
645 double rCut = 3.9*spacing;
646 double R=1.0;
647 Ghost<3, double> ghost(rCut);
648 // P V v_B RHS V_t P_anal RHS2 Polar cord
649 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);
650
651
652 auto &v_cl = create_vcluster();
653
654// openfpm::vector<aggregate<int>> bulk;
655// openfpm::vector<aggregate<int>> Surface;
656
657 auto it = Particles.getGridIterator(sz);
658 while (it.isNext()) {
659 auto key = it.get();
660 double x = -1.0+key.get(0) * it.getSpacing(0);
661 double y = -1.0+key.get(1) * it.getSpacing(1);
662 double z = -1.0+key.get(2) * it.getSpacing(2);
663 double r=sqrt(x*x+y*y+z*z);
664 if (r<R-spacing/2.0) {
665 Particles.add();
666 Particles.getLastPos()[0] = x;
667 Particles.getLastPos()[1] = y;
668 Particles.getLastPos()[2] = z;
669 Particles.getLastProp<8>()[0] = r;
670 if (r==0){
671 Particles.getLastProp<8>()[1] = 0.0;
672 }
673 else{
674 Particles.getLastProp<8>()[1] = std::atan2(sqrt(x*x+y*y),z);
675 }
676 Particles.getLastProp<8>()[2] = std::atan2(y,x);
677 }
678 ++it;
679 }
680
681 int n_sp=int(grd_sz)*int(grd_sz)*3;
682
683 double Golden_angle=M_PI * (3.0 - sqrt(5.0));
684
685 for(int i=1;i<=n_sp;i++)
686 {
687 double y = 1.0 - (i /double(n_sp - 1.0)) * 2.0;
688 double radius = sqrt(1 - y * y);
689 double Golden_theta = Golden_angle * i;
690 double x = cos(Golden_theta) * radius;
691 double z = sin(Golden_theta) * radius;
692
693 if (acos(z)==0 || acos(z)==M_PI){
694 std::cout<<"Theta 0/Pi "<<std::endl;
695 continue;
696 }
697
698 Particles.add();
699 Particles.getLastPos()[0] = x;
700 Particles.getLastPos()[1] = y;
701 Particles.getLastPos()[2] = z;
702 Particles.getLastProp<8>()[0] = 1.0 ;
703 Particles.getLastProp<8>()[1] = std::atan2(sqrt(x*x+y*y),z);
704 Particles.getLastProp<8>()[2] = std::atan2(y,x);
705 }
706 Particles.map();
707 Particles.ghost_get<0>();
708
709
710 std::unordered_map<const lm,double,key_hash,key_equal> Vr;
711 std::unordered_map<const lm,double,key_hash,key_equal> V1;
712 std::unordered_map<const lm,double,key_hash,key_equal> V2;
713 //Setting max mode l_max
714 constexpr int K = 2;
715 //Setting amplitudes to 0
716 for(int l=0;l<=K;l++){
717 for(int m=-l;m<=l;m++){
718 Vr[std::make_tuple(l,m)]=0.0;
719 V1[std::make_tuple(l,m)]=0.0;
720 V2[std::make_tuple(l,m)]=0.0;
721 }
722
723
724 }
725 //Setting some amplitude for boundary velocity
726 V1[std::make_tuple(1,0)]=1.0;
727
728 auto it2 = Particles.getDomainIterator();
729 while (it2.isNext()) {
730 auto p = it2.get();
731 Point<3, double> xp = Particles.getPos(p);
732 Point<3, double> xP = Particles.getProp<8>(p);
733 Particles.getProp<0>(p) =0;
734 if (xP[0]==1.0) {
735// Surface.add();
736// Surface.last().get<0>() = p.getKey();
737 Particles.getProp<0>(p) = 0;
738 std::vector<double> SVel;
739 SVel=openfpm::math::sumY<K>(xP[0],xP[1],xP[2],Vr,V1,V2);
740 double SP=openfpm::math::sumY_Scalar<K>(xP[0],xP[1],xP[2],Vr);
741 Particles.getProp<2>(p)[0] = SVel[0];
742 Particles.getProp<2>(p)[1] = SVel[1];
743 Particles.getProp<2>(p)[2] = SVel[2];
744 Particles.getProp<9>(p)[0] = SVel[0];
745 Particles.getProp<9>(p)[1] = SVel[1];
746 Particles.getProp<9>(p)[2] = SVel[2];
747 Particles.getProp<5>(p) = SP;
748 Particles.setSubset(p,1);
749
750 }
751 else {
752// bulk.add();
753// bulk.last().get<0>() = p.getKey();
754 Particles.setSubset(p,0);
755 Particles.getProp<0>(p) = 0;
756 Particles.getProp<1>(p)[0] = 0;
757 Particles.getProp<1>(p)[1] = 0;
758 Particles.getProp<1>(p)[2] = 0;
759 }
760 ++it2;
761 }
762
763 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);
764 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);
765
766 auto & bulk = Particles_bulk.getIds();
767 auto & Surface = Particles_surface.getIds();
768
769 for (int j = 0; j < bulk.size(); j++) {
770 auto p = bulk.get<0>(j);
771 Point<3, double> xp = Particles.getPos(p);
772 Point<3, double> xP = Particles.getProp<8>(p);
773
774 std::unordered_map<const lm,double,key_hash,key_equal> Ur;
775 std::unordered_map<const lm,double,key_hash,key_equal> U2;
776 std::unordered_map<const lm,double,key_hash,key_equal> U1;
777 std::unordered_map<const lm,double,key_hash,key_equal> Plm;
778
779 for (int l = 0; l <= K; l++) {
780 for (int m = -l; m <= l; m++) {
781 auto Er= Vr.find(std::make_tuple(l,m));
782 auto E1= V1.find(std::make_tuple(l,m));
783 auto E2= V2.find(std::make_tuple(l,m));
784 std::vector<double> Sol=openfpm::math::sph_anasol_u(nu,l,m,Er->second,E1->second,E2->second,xP[0]);
785 Ur[std::make_tuple(l,m)]=Sol[0];
786 U1[std::make_tuple(l,m)]=Sol[1];
787 U2[std::make_tuple(l,m)]=Sol[2];
788 Plm[std::make_tuple(l,m)]=Sol[3];
789 }
790
791 }
792
793 if(fabs(xP[0])>=1e-5 && xP[1]>1e-5 && (M_PI-xP[1])>=1e-5)
794 {
795 std::vector<double> SVel = openfpm::math::sumY<K>(xP[0], xP[1], xP[2], Ur, U1, U2);
796 Particles.getProp<9>(p)[0] = SVel[0];
797 Particles.getProp<9>(p)[1] = SVel[1];
798 Particles.getProp<9>(p)[2] = SVel[2];
799 Particles.getProp<5>(p) = openfpm::math::sumY_Scalar<K>(xP[0], xP[1], xP[2], Plm);
800 }
801 }
802
803 auto P = getV<0>(Particles);
804 auto V = getV<1>(Particles);
805 auto V_B = getV<2>(Particles);
806 V.setVarId(0);
807 auto DIV = getV<3>(Particles);
808 auto V_t = getV<4>(Particles);
809 auto P_anal = getV<5>(Particles);
810 auto temp=getV<6>(Particles);
811 auto RHS = getV<7>(Particles);
812 auto P_bulk = getV<0>(Particles_bulk);
813 auto RHS_bulk = getV<7>(Particles_bulk);
814 auto V_anal = getV<9>(Particles);
815
816 V_t=V;
817 P=0;
818 P_bulk=0;
819 eq_id vx,vy,vz;
820
821 vx.setId(0);
822 vy.setId(1);
823 vz.setId(2);
824
825 double sampling=3.1;
826 double sampling2=1.9;
827 double rCut2=3.9*spacing;
828
829 Derivative_x Dx(Particles, 2, rCut,sampling, support_options::RADIUS),B_Dx(Particles_bulk, 2, rCut,sampling, support_options::RADIUS);
830 Derivative_y Dy(Particles, 2, rCut,sampling, support_options::RADIUS),B_Dy(Particles_bulk, 2, rCut,sampling, support_options::RADIUS);
831 Derivative_z Dz(Particles, 2, rCut,sampling, support_options::RADIUS),B_Dz(Particles_bulk, 2, rCut,sampling, support_options::RADIUS);
832 Derivative_xx Dxx(Particles, 2, rCut2,sampling2,support_options::RADIUS);
833 Derivative_yy Dyy(Particles, 2, rCut2,sampling2,support_options::RADIUS);
834 Derivative_zz Dzz(Particles, 2, rCut2,sampling2,support_options::RADIUS);
835
836 //std::cout << "DCPSE KERNELS DONE" << std::endl;
837 petsc_solver<double> solverPetsc;
838 solverPetsc.setPreconditioner(PCNONE);
839 timer tt;
840 double sum=0,sum1=0;
841 V_t=V;
842 double V_err_eps = 1e-5;
843
844 double V_err = 1, V_err_old;
845 int n = 0;
846 int nmax = 30;
847 int ctr = 0, errctr, Vreset = 0;
848 V_err = 1;
849 n = 0;
850 double solvetime=0;
851 while (V_err >= V_err_eps && n <= nmax) {
852 //Particles.write_frame("StokesSphere",n);
853 Particles.ghost_get<0>(SKIP_LABELLING);
854 RHS_bulk[0] = B_Dx(P);
855 RHS_bulk[1] = B_Dy(P);
856 RHS_bulk[2] = B_Dz(P);
857 DCPSE_scheme<equations3d3, decltype(Particles)> Solver(Particles);
858 auto Stokes1 = nu * (Dxx(V[0])+Dyy(V[0])+Dzz(V[0]));
859 auto Stokes2 = nu * (Dxx(V[1])+Dyy(V[1])+Dzz(V[1]));
860 auto Stokes3 = nu * (Dxx(V[2])+Dyy(V[2])+Dzz(V[2]));
861 Solver.impose(Stokes1, bulk, RHS[0], vx);
862 Solver.impose(Stokes2, bulk, RHS[1], vy);
863 Solver.impose(Stokes3, bulk, RHS[2], vz);
864 Solver.impose(V[0], Surface, V_B[0], vx);
865 Solver.impose(V[1], Surface, V_B[1], vy);
866 Solver.impose(V[2], Surface, V_B[2], vz);
867 Solver.impose_x_ig(bulk, V[0], vx);
868 Solver.impose_x_ig(bulk, V[1], vy);
869 Solver.impose_x_ig(bulk, V[2], vz);
870 Solver.impose_x_ig(Surface, V[0], vx);
871 Solver.impose_x_ig(Surface, V[1], vy);
872 Solver.impose_x_ig(Surface, V[2], vz);
873 tt.start();
874 Solver.solve_with_solver_ig(solverPetsc, V[0], V[1], V[2]);
875 tt.stop();
876 solvetime+=tt.getwct();
877 //Solver.solve(V[0],V[1],V[2]);
878 //std::cout << "Stokes Solved" << std::endl;
879 Particles.ghost_get<1>();
880 DIV = -(Dx(V[0])+Dy(V[1])+Dz(V[2]));
881 P_bulk = P + DIV;
882 sum = 0;
883 sum1 = 0;
884 for (int j = 0; j < bulk.size(); j++) {
885 auto p = bulk.get<0>(j);
886 sum += (Particles.getProp<4>(p)[0] - Particles.getProp<1>(p)[0]) *
887 (Particles.getProp<4>(p)[0] - Particles.getProp<1>(p)[0]) +
888 (Particles.getProp<4>(p)[1] - Particles.getProp<1>(p)[1]) *
889 (Particles.getProp<4>(p)[1] - Particles.getProp<1>(p)[1]) +
890 (Particles.getProp<4>(p)[2] - Particles.getProp<1>(p)[2]) *
891 (Particles.getProp<4>(p)[2] - Particles.getProp<1>(p)[2]);
892 sum1 += Particles.getProp<1>(p)[0] * Particles.getProp<1>(p)[0] +
893 Particles.getProp<1>(p)[1] * Particles.getProp<1>(p)[1] +
894 Particles.getProp<1>(p)[2] * Particles.getProp<1>(p)[2];
895 }
896 sum = sqrt(sum);
897 sum1 = sqrt(sum1);
898 v_cl.sum(sum);
899 v_cl.sum(sum1);
900 v_cl.execute();
901 V_t = V;
902 Particles.ghost_get<1>(SKIP_LABELLING);
903 V_err_old = V_err;
904 V_err = sum / sum1;
905 if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-14) {
906 errctr++;
907 } else {
908 errctr = 0;
909 }
910 if (n > 3) {
911 if (errctr > 1) {
912 std::cout << "CONVERGENCE LOOP BROKEN DUE TO INCREASE/VERY SLOW DECREASE IN ERROR" << std::endl;
913 Vreset = 1;
914 break;
915 } else {
916 Vreset = 0;
917 }
918 }
919 n++;
920
921 }
922 //std::cout << "Total Solver time (wct):"<<solvetime<< std::endl;
923
924 V_t=0;
925
926 double worst=0;
927 double L2=0;
928 for (int j = 0; j < bulk.size(); j++) {
929 auto p = bulk.get<0>(j);
930 Point<3,double> xP=Particles.getProp<8>(p);
931 if(xP[0]>=1e-5 && xP[1]>1e-5 && (M_PI-xP[1])>=1e-5)
932 {
933 double dx=Particles.getProp<9>(p)[0] - Particles.getProp<1>(p)[0];
934 double dy=Particles.getProp<9>(p)[1] - Particles.getProp<1>(p)[1];
935 double dz=Particles.getProp<9>(p)[2] - Particles.getProp<1>(p)[2];
936 Particles.getProp<4>(p)[0]=fabs(dx);
937 Particles.getProp<4>(p)[1]=fabs(dy);
938 Particles.getProp<4>(p)[2]=fabs(dz);
939 L2 += dx*dx+dy*dy+dz*dz;
940 if (std::max({fabs(dx),fabs(dy),fabs(dz)}) > worst) {
941 worst = std::max({fabs(dx),fabs(dy),fabs(dz)});
942 }
943
944 }
945 }
946
947 v_cl.sum(worst);
948 v_cl.sum(L2);
949 v_cl.execute();
950
951 /* if (v_cl.rank() == 0) {
952 std::cout<<"Gd,Surf,Bulk Size: "<<grd_sz<<","<<Surface.size()<<","<<bulk.size()<<std::endl;
953 std::cout << "L2_Final: " <<sqrt(L2)<<","<<sqrt(L2/(bulk.size()+Surface.size()))
954 << std::endl;
955 std::cout << "L_inf_Final: " << worst
956 << std::endl;
957 }
958 std::cout << "L_inf_Final_test: " << worst;*/
959 //Particles.write("StokesSphere");
960 BOOST_REQUIRE(worst<1e-3);
961
962 }
963BOOST_AUTO_TEST_SUITE_END()
964
965#endif
966#endif
Box
This class represent an N-dimensional box.
Definition Box.hpp:61
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:203
openfpm::vector::size
size_t size()
Stub size.
Definition map_vector.hpp:211
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::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:305
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