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