OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
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 
29 BOOST_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 
615 BOOST_AUTO_TEST_SUITE_END()
616 #endif
617 #endif
vector_dist_subset
Definition: vector_dist_subset.hpp:69
vector_dist_ws
Definition: vector_dist_subset.hpp:16
Point
This class implement the point shape in an N-dimensional space.
Definition: Point.hpp:27
openfpm::vector::size
size_t size()
Stub size.
Definition: map_vector.hpp:211
Ghost
Definition: Ghost.hpp:39
petsc_solver< double >
This class is able to do Matrix inversion in parallel with PETSC solvers.
Definition: petsc_solver.hpp:101
Box
This class represent an N-dimensional box.
Definition: Box.hpp:60
vector_dist
Distributed vector.
Definition: vector_dist.hpp:297
aggregate
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Definition: aggregate.hpp:214
sum
It model an expression expr1 + ... exprn.
Definition: sum.hpp:92
Point_test
Test structure used for several test.
Definition: Point_test.hpp:105
openfpm::vector
Implementation of 1-D std::vector like structure.
Definition: map_vector.hpp:202
eq_id
Definition: eq_solve_common.hpp:14