OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
DCPSE_op_subset_test.cpp
1/*
2 * DCPSE_op_test.cpp
3 *
4 * Created on: May 15, 2020
5 * Author: Abhinav Singh
6 *
7 */
8#include "config.h"
9#ifdef HAVE_EIGEN
10#ifdef HAVE_PETSC
11
12
13#define BOOST_TEST_DYN_LINK
14
15#include "util/util_debug.hpp"
16#include <boost/test/unit_test.hpp>
17#include <iostream>
18#include "../DCPSE_op.hpp"
19#include "../DCPSE_Solver.hpp"
20#include "Operators/Vector/vector_dist_operators.hpp"
21#include "Vector/vector_dist_subset.hpp"
22#include "../EqnsStruct.hpp"
23
24//template<typename T>
25//struct Debug;
26
27
28
29BOOST_AUTO_TEST_SUITE(dcpse_op_subset_suite_tests)
30
31 BOOST_AUTO_TEST_CASE(dcpse_op_subset_tests) {
32// int rank;
33// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
34 size_t edgeSemiSize = 40;
35 const size_t sz[2] = {2 * edgeSemiSize, 2 * edgeSemiSize};
36 Box<2, double> box({0, 0}, {1.0,1.0});
37 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
38 double spacing[2];
39 spacing[0] = 1.0 / (sz[0] - 1);
40 spacing[1] = 1.0 / (sz[1] - 1);
41 double rCut = 3.9 * spacing[0];
42 int ord = 2;
43 double sampling_factor = 4.0;
44 Ghost<2, double> ghost(rCut);
45 BOOST_TEST_MESSAGE("Init vector_dist...");
46 double sigma2 = spacing[0] * spacing[1] / (2 * 4);
47
48 vector_dist_ws<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double>> Particles(0, box,
49 bc,
50 ghost);
51 Particles.setPropNames({"0","1","2","3","4","5","6"});
52 //Init_DCPSE(Particles)
53 BOOST_TEST_MESSAGE("Init Particles...");
54 std::mt19937 rng{6666666};
55
56 std::normal_distribution<> gaussian{0, sigma2};
57
58// openfpm::vector<aggregate<int>> bulk;
59// openfpm::vector<aggregate<int>> boundary;
60
61 auto it = Particles.getGridIterator(sz);
62 size_t pointId = 0;
63 size_t counter = 0;
64 double minNormOne = 999;
65 while (it.isNext())
66 {
67 Particles.add();
68 auto key = it.get();
69 mem_id k0 = key.get(0);
70 double x = k0 * spacing[0];
71 Particles.getLastPos()[0] = x;//+ gaussian(rng);
72 mem_id k1 = key.get(1);
73 double y = k1 * spacing[1];
74 Particles.getLastPos()[1] = y;//+gaussian(rng);
75 // Here fill the function value
76 Particles.template getLastProp<0>() = sin(Particles.getLastPos()[0]) + sin(Particles.getLastPos()[1]);
77
78 if (k0 != 0 && k1 != 0 && k0 != sz[0] -1 && k1 != sz[1] - 1)
79 {
80// bulk.add();
81// bulk.template get<0>(bulk.size()-1) = Particles.size_local() - 1;
82
83 Particles.getLastSubset(0);
84 }
85 else
86 {
87// boundary.add();
88// boundary.template get<0>(boundary.size()-1) = Particles.size_local() - 1;
89 Particles.getLastSubset(1);
90 }
91
92
93 ++counter;
94 ++it;
95 }
96 BOOST_TEST_MESSAGE("Sync Particles across processors...");
97
98 Particles.map();
99 Particles.ghost_get<0>();
100
101 auto git = Particles.getGhostIterator();
102
103 while (git.isNext())
104 {
105 auto p = git.get();
106
107 Particles.template getProp<0>(p) = std::numeric_limits<double>::quiet_NaN();
108
109 ++git;
110 }
111
112 vector_dist_subset<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double>> Particles_bulk(Particles,0);
113 vector_dist_subset<2, double, aggregate<double, double, double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double>> Particles_boundary(Particles,1);
114 auto & boundary = Particles_boundary.getIds();
115
116 // move particles
117 auto P = getV<0>(Particles);
118 auto Out = getV<1>(Particles);
119 auto Pb = getV<2>(Particles);
120 auto Out_V = getV<3>(Particles);
121
122
123 auto P_bulk = getV<2>(Particles_bulk);
124 auto Out_bulk = getV<1>(Particles_bulk);
125 auto Out_V_bulk = getV<3>(Particles_bulk);
126 Out=10;
127 P_bulk = 5;
128
129 P_bulk=Pb+Out;
130// Particles.write("Test_output_subset");
131
132 // Create the subset
133
134 /* Derivative_x Dx(Particles, 2, rCut);
135 Derivative_y Dy(Particles, 2, rCut);
136 Derivative_x Dx_bulk(Particles_bulk, 2, rCut);
137*/
138 Derivative_x Dx_bulk(Particles_bulk, 2, rCut,sampling_factor, support_options::RADIUS);
139 Derivative_y Dy_bulk(Particles_bulk, 2, rCut,sampling_factor, support_options::RADIUS);
140
141 Out_bulk = Dx_bulk(P);
142 Out_V_bulk[0] = P + Dx_bulk(P);
143 Out_V_bulk[1] = Out_V[0] +Dy_bulk(P);
144
145 // Check
146 bool is_nan = false;
147
148 auto & v_cl = create_vcluster();
149 if (v_cl.size() > 1)
150 {
151 auto it2 = Particles_bulk.getDomainIterator();
152 while (it2.isNext())
153 {
154 auto p = it2.get();
155
156 /* BOOST_REQUIRE_EQUAL(Particles_bulk.getProp<2>(p),15.0);
157 BOOST_REQUIRE(fabs(Particles_bulk.getProp<1>(p) - cos(Particles_bulk.getPos(p)[0])) < 0.005 );
158 BOOST_REQUIRE(fabs(Particles_bulk.getProp<3>(p)[0] - Particles_bulk.getProp<0>(p) - cos(Particles_bulk.getPos(p)[0])) < 0.001 );
159 BOOST_REQUIRE(fabs(Particles_bulk.getProp<3>(p)[1] - Particles_bulk.getProp<3>(p)[0] - cos(Particles_bulk.getPos(p)[1])) < 0.001 );*/
160
161 is_nan |= std::isnan(Particles_bulk.template getProp<1>(p));
162
163 // Particles_bulk.template getProp<0>(p) = fabs(Particles_bulk.getProp<1>(p) - cos(Particles_bulk.getPos(p)[0]));
164
165 ++it2;
166 }
167
168 BOOST_REQUIRE_EQUAL(is_nan,true);
169 }
170
171// P_bulk = Dx_bulk(P_bulk); <------------ Incorrect produce error message
172// P = Dx_bulk(P); <------- Incorrect produce overflow
173
174 Particles.ghost_get<0>();
175 Particles.write("TEST");
176
177 for (int i = 0 ; i < boundary.size() ; i++)
178 {
179 Particles.template getProp<0>(boundary.template get<0>(i)) = std::numeric_limits<double>::quiet_NaN();
180 }
181
182 Particles.ghost_get<0>();
183
184 Out_bulk = Dx_bulk(P);
185 Out_V_bulk[0] = P + Dx_bulk(P);
186 Out_V_bulk[1] = Out_V[0] +Dy_bulk(P);
187
188 auto it2 = Particles_bulk.getDomainIterator();
189 while (it2.isNext())
190 {
191 auto p = it2.get();
192
193 BOOST_REQUIRE_EQUAL(Particles_bulk.getProp<2>(p),15.0);
194 BOOST_REQUIRE(fabs(Particles_bulk.getProp<1>(p) - cos(Particles_bulk.getPos(p)[0])) < 0.005 );
195 BOOST_REQUIRE(fabs(Particles_bulk.getProp<3>(p)[0] - Particles_bulk.getProp<0>(p) - cos(Particles_bulk.getPos(p)[0])) < 0.001 );
196 BOOST_REQUIRE(fabs(Particles_bulk.getProp<3>(p)[1] - Particles_bulk.getProp<3>(p)[0] - cos(Particles_bulk.getPos(p)[1])) < 0.001 );
197
198 ++it2;
199 }
200
201
202 }
203
204 BOOST_AUTO_TEST_CASE(dcpse_op_subset_PC_lid) {
205// int rank;
206// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
207 constexpr int x = 0;
208 constexpr int y = 1;
209 size_t edgeSemiSize = 20;
210 const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
211 Box<2, double> box({0, 0}, {1.0,1.0});
212 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
213 double spacing[2];
214 spacing[0] = 1.0 / (sz[0] - 1);
215 spacing[1] = 1.0 / (sz[1] - 1);
216 double rCut = 3.9 * spacing[0];
217 int ord = 2;
218 double sampling_factor = 4.0;
219 Ghost<2, double> ghost(rCut);
220 BOOST_TEST_MESSAGE("Init vector_dist...");
221 double sigma2 = spacing[0] * spacing[1] / (2 * 4);
222 auto &v_cl = create_vcluster();
223
224 typedef aggregate<double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,VectorS<2, double>,double> particle_type;
225
226 vector_dist_ws<2, double, particle_type> Particles(0, box,bc,ghost);
227
228 //Init_DCPSE(Particles)
229 BOOST_TEST_MESSAGE("Init Particles...");
230
231// openfpm::vector<aggregate<int>> bulk;
232// openfpm::vector<aggregate<int>> boundary;
233
234 auto it = Particles.getGridIterator(sz);
235 while (it.isNext())
236 {
237 Particles.add();
238 auto key = it.get();
239 mem_id k0 = key.get(0);
240 double xp0 = k0 * spacing[0];
241 Particles.getLastPos()[0] = xp0;
242 mem_id k1 = key.get(1);
243 double yp0 = k1 * spacing[1];
244 Particles.getLastPos()[1] = yp0;
245 ++it;
246 }
247 BOOST_TEST_MESSAGE("Sync Particles across processors...");
248 Particles.map();
249 Particles.ghost_get<0>();
250
251 auto it2 = Particles.getDomainIterator();
252 while (it2.isNext()) {
253 auto p = it2.get();
254 Point<2, double> xp = Particles.getPos(p);
255 if (xp[0] != 0 && xp[1] != 0 && xp[0] != 1.0 && xp[1] != 1.0) {
256// bulk.add();
257// bulk.last().get<0>() = p.getKey();
258 Particles.setSubset(p,0);
259 Particles.getProp<3>(p)[x] = 3.0;
260 Particles.getProp<3>(p)[y] = 3.0;
261 } else {
262// boundary.add();
263// boundary.last().get<0>() = p.getKey();
264 Particles.setSubset(p,1);
265 Particles.getProp<3>(p)[x] = xp[0]*xp[0]+xp[1]*xp[1];
266 Particles.getProp<3>(p)[y] = xp[0]*xp[0]-2*xp[0]*xp[1];
267 }
268 Particles.getProp<6>(p)[x] = xp[0]*xp[0]+xp[1]*xp[1];
269 Particles.getProp<6>(p)[y] = xp[0]*xp[0]-2*xp[0]*xp[1];
270 Particles.getProp<7>(p) = xp[0]+xp[1]-1.0;
271
272 ++it2;
273 }
274
275 vector_dist_subset<2, double, particle_type> Particles_bulk(Particles,0);
276 vector_dist_subset<2, double, particle_type> Particles_boundary(Particles,1);
277 auto & bulk = Particles_bulk.getIds();
278 auto & boundary = Particles_boundary.getIds();
279
280 auto P = getV<0>(Particles);
281 auto V = getV<1>(Particles);
282 auto RHS = getV<2>(Particles);
283 auto dV = getV<3>(Particles);
284 auto div = getV<4>(Particles);
285 auto V_star = getV<5>(Particles);
286
287
288 auto P_bulk = getV<0>(Particles_bulk);
289 auto RHS_bulk =getV<2>(Particles_bulk);
290
291 P_bulk = 0;
292
293 Derivative_x Dx(Particles, 2, rCut,sampling_factor, support_options::RADIUS);
294 Derivative_xx Dxx(Particles, 2, rCut,sampling_factor, support_options::RADIUS);
295 Derivative_yy Dyy(Particles, 2, rCut,sampling_factor, support_options::RADIUS);
296 Derivative_y Dy(Particles, 2, rCut,sampling_factor, support_options::RADIUS);
297 Derivative_x Bulk_Dx(Particles_bulk, 2, rCut,sampling_factor, support_options::RADIUS);
298 Derivative_y Bulk_Dy(Particles_bulk, 2, rCut,sampling_factor, support_options::RADIUS);
299
300 int n = 0, nmax = 5, ctr = 0, errctr=1, Vreset = 0;
301 double V_err=1;
302 if (Vreset == 1) {
303 P_bulk = 0;
304 P = 0;
305 Vreset = 0;
306 }
307 P=0;
308 eq_id vx,vy;
309 vx.setId(0);
310 vy.setId(1);
311 double sum, sum1, sum_k,V_err_eps=1e-3,V_err_old;
312 auto Stokes1=Dxx(V[x])+Dyy(V[x]);
313 auto Stokes2=Dxx(V[y])+Dyy(V[y]);
314
315 petsc_solver<double> solverPetsc;
316 //solverPetsc.setSolver(KSPGMRES);
317 //solverPetsc.setRestart(250);
318 //solverPetsc.setPreconditioner(PCJACOBI);
319 V_star=0;
320 RHS[x] = dV[x];
321 RHS[y] = dV[y];
322 while (V_err >= V_err_eps && n <= nmax) {
323 Particles.ghost_get<0>(SKIP_LABELLING);
324 RHS_bulk[x] = dV[x] + Bulk_Dx(P);
325 RHS_bulk[y] = dV[y] + Bulk_Dy(P);
326 DCPSE_scheme<equations2d2, decltype(Particles)> Solver(Particles);
327 Solver.impose(Stokes1, bulk, RHS[0], vx);
328 Solver.impose(Stokes2, bulk, RHS[1], vy);
329 Solver.impose(V[x], boundary, RHS[0], vx);
330 Solver.impose(V[y], boundary, RHS[1], vy);
331
332 /*auto A=Solver.getA(options_solver::STANDARD);
333 //A.getMatrixTriplets().save("Tripletes");
334 A.write("Mat_lid");*/
335
336 Solver.solve_with_solver(solverPetsc, V[x], V[y]);
337 Particles.ghost_get<1>(SKIP_LABELLING);
338 div = -(Dx(V[x]) + Dy(V[y]));
339 P_bulk = P + div;
340 sum = 0;
341 sum1 = 0;
342
343 for (int j = 0; j < bulk.size(); j++) {
344 auto p = bulk.get<0>(j);
345 sum += (Particles.getProp<5>(p)[0] - Particles.getProp<1>(p)[0]) *
346 (Particles.getProp<5>(p)[0] - Particles.getProp<1>(p)[0]) +
347 (Particles.getProp<5>(p)[1] - Particles.getProp<1>(p)[1]) *
348 (Particles.getProp<5>(p)[1] - Particles.getProp<1>(p)[1]);
349 sum1 += Particles.getProp<1>(p)[0] * Particles.getProp<1>(p)[0] +
350 Particles.getProp<1>(p)[1] * Particles.getProp<1>(p)[1];
351 }
352
353 sum = sqrt(sum);
354 sum1 = sqrt(sum1);
355 V_star = V;
356 v_cl.sum(sum);
357 v_cl.sum(sum1);
358 v_cl.execute();
359 V_err_old = V_err;
360 V_err = sum / sum1;
361 if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-8) {
362 errctr++;
363 //alpha_P -= 0.1;
364 } else {
365 errctr = 0;
366 }
367 if (n > 3) {
368 if (errctr > 3) {
369 std::cout << "CONVERGENCE LOOP BROKEN DUE TO INCREASE/VERY SLOW DECREASE IN ERROR" << std::endl;
370 Vreset = 1;
371 break;
372 } else {
373 Vreset = 0;
374 }
375 }
376 n++;
377 if (v_cl.rank() == 0) {
378 std::cout << "Rel l2 cgs err in V = " << V_err << " at " << n << std::endl;
379 }
380 }
381 double worst1 = 0.0;
382 double worst2 = 0.0;
383
384 for(int j=0;j<bulk.size();j++)
385 { auto p=bulk.get<0>(j);
386 if (fabs(Particles.getProp<6>(p)[0] - Particles.getProp<1>(p)[0]) >= worst1) {
387 worst1 = fabs(Particles.getProp<6>(p)[0] - Particles.getProp<1>(p)[0]);
388 }
389 }
390 for(int j=0;j<bulk.size();j++)
391 { auto p=bulk.get<0>(j);
392 if (fabs(Particles.getProp<6>(p)[1] - Particles.getProp<1>(p)[1]) >= worst2) {
393 worst2 = fabs(Particles.getProp<6>(p)[1] - Particles.getProp<1>(p)[1]);
394 }
395 }
396 //Particles.deleteGhost();
397 //Particles.write("PC_subset_lid");
398 std::cout << "Maximum Analytic Error in Vx: " << worst1 << std::endl;
399 std::cout << "Maximum Analytic Error in Vy: " << worst2 << std::endl;
400 BOOST_REQUIRE(worst1 < 0.03);
401 BOOST_REQUIRE(worst2 < 0.03);
402
403 }
404
405
406 BOOST_AUTO_TEST_CASE(dcpse_op_subset_PC_lid2) {
407// int rank;
408// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
409 constexpr int x = 0;
410 constexpr int y = 1;
411 size_t edgeSemiSize = 20;
412 const size_t sz[2] = {2 * edgeSemiSize+1, 2 * edgeSemiSize+1};
413 Box<2, double> box({0, 0}, {1.0,1.0});
414 size_t bc[2] = {NON_PERIODIC, NON_PERIODIC};
415 double spacing[2];
416 spacing[0] = 1.0 / (sz[0] - 1);
417 spacing[1] = 1.0 / (sz[1] - 1);
418 double rCut = 3.9 * spacing[0];
419 int ord = 2;
420 double sampling_factor = 4.0;
421 Ghost<2, double> ghost(rCut);
422 BOOST_TEST_MESSAGE("Init vector_dist...");
423 auto &v_cl = create_vcluster();
424
425 vector_dist<2, double, aggregate<double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,VectorS<2, double>,double>> Particles(0, box,
426 bc,
427 ghost);
428 vector_dist<2, double, aggregate<double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,VectorS<2, double>,double>> Particles_subset(Particles.getDecomposition(), 0);
429
430
431 //Init_DCPSE(Particles)
432 BOOST_TEST_MESSAGE("Init Particles...");
433
434 openfpm::vector<aggregate<int>> bulk;
435 openfpm::vector<aggregate<int>> boundary;
436
437 auto it = Particles.getGridIterator(sz);
438 size_t pointId = 0;
439 double minNormOne = 999;
440 while (it.isNext())
441 {
442 Particles.add();
443 auto key = it.get();
444 mem_id k0 = key.get(0);
445 double xp0 = k0 * spacing[0];
446 Particles.getLastPos()[0] = xp0;
447 mem_id k1 = key.get(1);
448 double yp0 = k1 * spacing[1];
449 Particles.getLastPos()[1] = yp0;
450 ++it;
451 }
452 BOOST_TEST_MESSAGE("Sync Particles across processors...");
453 Particles.map();
454 Particles.ghost_get<0>();
455 auto it2 = Particles.getDomainIterator();
456 while (it2.isNext()) {
457 auto p = it2.get();
458 Point<2, double> xp = Particles.getPos(p);
459 if (xp[0] != 0 && xp[1] != 0 && xp[0] != 1.0 && xp[1] != 1.0) {
460 bulk.add();
461 bulk.last().get<0>() = p.getKey();
462 Particles.getProp<3>(p)[x] = 3.0;
463 Particles.getProp<3>(p)[y] = 3.0;
464 } else {
465 boundary.add();
466 boundary.last().get<0>() = p.getKey();
467 Particles.getProp<3>(p)[x] = xp[0]*xp[0]+xp[1]*xp[1];
468 Particles.getProp<3>(p)[y] = xp[0]*xp[0]-2*xp[0]*xp[1];
469 }
470 Particles.getProp<6>(p)[x] = xp[0]*xp[0]+xp[1]*xp[1];
471 Particles.getProp<6>(p)[y] = xp[0]*xp[0]-2*xp[0]*xp[1];
472 Particles.getProp<7>(p) = xp[0]+xp[1]-1.0;
473
474 ++it2;
475 }
476
477 for (int i = 0; i < bulk.size(); i++) {
478 Particles_subset.add();
479 Particles_subset.getLastPos()[0] = Particles.getPos(bulk.template get<0>(i))[0];
480 Particles_subset.getLastPos()[1] = Particles.getPos(bulk.template get<0>(i))[1];
481 }
482 Particles_subset.map();
483 Particles_subset.ghost_get<0>();
484
485
486
487 auto P = getV<0>(Particles);
488 auto V = getV<1>(Particles);
489 auto RHS = getV<2>(Particles);
490 auto dV = getV<3>(Particles);
491 auto div = getV<4>(Particles);
492 auto V_star = getV<5>(Particles);
493
494 auto P_bulk = getV<0>(Particles_subset);
495 auto Grad_bulk= getV<2>(Particles_subset);
496
497 P_bulk = 0;
498
499 Derivative_x Dx(Particles, 2, rCut,sampling_factor, support_options::RADIUS);
500 Derivative_x Bulk_Dx(Particles_subset, 2, rCut,sampling_factor, support_options::RADIUS);
501 Derivative_xx Dxx(Particles, 2, rCut,sampling_factor, support_options::RADIUS);
502 Derivative_yy Dyy(Particles, 2, rCut,sampling_factor, support_options::RADIUS);
503 Derivative_y Dy(Particles, 2, rCut,sampling_factor, support_options::RADIUS),Bulk_Dy(Particles_subset, 2, rCut,sampling_factor, support_options::RADIUS);;
504
505 int n = 0, nmax = 5, ctr = 0, errctr=0, Vreset = 0;
506 double V_err=1;
507 if (Vreset == 1) {
508 P_bulk = 0;
509 P = 0;
510 Vreset = 0;
511 }
512 P=0;
513 eq_id vx,vy;
514 vx.setId(0);
515 vy.setId(1);
516 double sum, sum1, sum_k,V_err_eps=1e-3,V_err_old;
517 auto Stokes1=Dxx(V[x])+Dyy(V[x]);
518 auto Stokes2=Dxx(V[y])+Dyy(V[y]);
519
520 petsc_solver<double> solverPetsc;
521 //solverPetsc.setSolver(KSPGMRES);
522 //solverPetsc.setRestart(250);
523 //solverPetsc.setPreconditioner(PCJACOBI);
524 V_star=0;
525 while (V_err >= V_err_eps && n <= nmax) {
526 RHS[x] = dV[x];
527 RHS[y] = dV[y];
528 Particles_subset.ghost_get<0>(SKIP_LABELLING);
529
530 Grad_bulk[x] = Bulk_Dx(P_bulk);
531 Grad_bulk[y] = Bulk_Dy(P_bulk);
532 for (int i = 0; i < bulk.size(); i++) {
533 Particles.template getProp<2>(bulk.template get<0>(i))[x] += Particles_subset.getProp<2>(i)[x];
534 Particles.template getProp<2>(bulk.template get<0>(i))[y] += Particles_subset.getProp<2>(i)[y];
535 }
536
537 DCPSE_scheme<equations2d2, decltype(Particles)> Solver(Particles);
538 Solver.impose(Stokes1, bulk, RHS[0], vx);
539 Solver.impose(Stokes2, bulk, RHS[1], vy);
540 Solver.impose(V[x], boundary, RHS[0], vx);
541 Solver.impose(V[y], boundary, RHS[1], vy);
542 Solver.solve_with_solver(solverPetsc, V[x], V[y]);
543 Particles.ghost_get<1>(SKIP_LABELLING);
544 div = -(Dx(V[x]) + Dy(V[y]));
545 P = P + div;
546 for (int i = 0; i < bulk.size(); i++) {
547 Particles_subset.getProp<0>(i) = Particles.template getProp<0>(bulk.template get<0>(i));
548 }
549 sum = 0;
550 sum1 = 0;
551
552 for (int j = 0; j < bulk.size(); j++) {
553 auto p = bulk.get<0>(j);
554 sum += (Particles.getProp<5>(p)[0] - Particles.getProp<1>(p)[0]) *
555 (Particles.getProp<5>(p)[0] - Particles.getProp<1>(p)[0]) +
556 (Particles.getProp<5>(p)[1] - Particles.getProp<1>(p)[1]) *
557 (Particles.getProp<5>(p)[1] - Particles.getProp<1>(p)[1]);
558 sum1 += Particles.getProp<1>(p)[0] * Particles.getProp<1>(p)[0] +
559 Particles.getProp<1>(p)[1] * Particles.getProp<1>(p)[1];
560 }
561
562 sum = sqrt(sum);
563 sum1 = sqrt(sum1);
564 V_star=V;
565 v_cl.sum(sum);
566 v_cl.sum(sum1);
567 v_cl.execute();
568 V_err_old = V_err;
569 V_err = sum / sum1;
570 if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-8) {
571 errctr++;
572 //alpha_P -= 0.1;
573 } else {
574 errctr = 0;
575 }
576 if (n > 3) {
577 if (errctr > 3) {
578 std::cout << "CONVERGENCE LOOP BROKEN DUE TO INCREASE/VERY SLOW DECREASE IN ERROR" << std::endl;
579 Vreset = 1;
580 break;
581 } else {
582 Vreset = 0;
583 }
584 }
585 n++;
586 if (v_cl.rank() == 0) {
587 std::cout << "Rel l2 cgs err in V = " << V_err << " at " << n << std::endl;
588
589 }
590 }
591 double worst1 = 0.0;
592 double worst2 = 0.0;
593
594 for(int j=0;j<bulk.size();j++)
595 { auto p=bulk.get<0>(j);
596 if (fabs(Particles.getProp<6>(p)[0] - Particles.getProp<1>(p)[0]) >= worst1) {
597 worst1 = fabs(Particles.getProp<6>(p)[0] - Particles.getProp<1>(p)[0]);
598 }
599 }
600 for(int j=0;j<bulk.size();j++)
601 { auto p=bulk.get<0>(j);
602 if (fabs(Particles.getProp<6>(p)[1] - Particles.getProp<1>(p)[1]) >= worst2) {
603 worst2 = fabs(Particles.getProp<6>(p)[1] - Particles.getProp<1>(p)[1]);
604 }
605 }
606
607 std::cout << "Maximum Analytic Error in slice x: " << worst1 << std::endl;
608 std::cout << "Maximum Analytic Error in slice y: " << worst2 << std::endl;
609 BOOST_REQUIRE(worst1 < 0.03);
610 BOOST_REQUIRE(worst2 < 0.03);
611
612 //Particles.write("PC_subset_lid2");
613 }
614
615BOOST_AUTO_TEST_SUITE_END()
616#endif
617#endif
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
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
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
aggregate
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Definition aggregate.hpp:215
sum
It model an expression expr1 + ... exprn.
Definition sum.hpp:93