OpenFPM  5.2.0
Project that contain the implementation of distributed structures
FD_Solver_test.cpp
1 #include "config.h"
2 #ifdef HAVE_EIGEN
3 #ifdef HAVE_PETSC
4 
5 #define BOOST_TEST_DYN_LINK
6 #include <boost/test/unit_test.hpp>
7 
8 #include "FD_Solver.hpp"
9 #include "Solvers/petsc_solver.hpp"
10 #include "Solvers/umfpack_solver.hpp"
11 #include "FD_expressions.hpp"
12 #include "FD_op.hpp"
13 #include "Grid/staggered_dist_grid.hpp"
14 #include "util/EqnsStructFD.hpp"
15 
16 
17 BOOST_AUTO_TEST_SUITE( FD_Solver_test )
18 
19 
20 BOOST_AUTO_TEST_CASE(solver_check_diagonal)
21 {
22  const size_t sz[2] = {81,81};
23  Box<2, double> box({0, 0}, {1, 1});
24  periodicity<2> bc = {NON_PERIODIC, NON_PERIODIC};
25  Ghost<2,long int> ghost(1);
26 
27  grid_dist_id<2, double, aggregate<double,double,double>> domain(sz, box, ghost, bc);
28 
29 
30  auto it = domain.getDomainIterator();
31  while (it.isNext())
32  {
33  auto key = it.get();
34  auto gkey = it.getGKey(key);
35  double x = gkey.get(0) * domain.spacing(0);
36  double y = gkey.get(1) * domain.spacing(1);
37  domain.get<1>(key) = x+y;
38 
39  ++it;
40  }
41 
42  domain.ghost_get<0>();
43  petsc_solver<double> pet_sol;
44  pet_sol.setPreconditioner(PCNONE);
45 
46  auto v = FD::getV<0>(domain);
47  auto RHS= FD::getV<1>(domain);
48  auto sol= FD::getV<2>(domain);
49  FD::LInfError LInfError;
50  FD_scheme<equations2d1,decltype(domain)> Solver(ghost,domain);
51  Solver.impose(5.0*v,{0,0},{80,80}, prop_id<1>());
52  Solver.solve_with_solver(pet_sol,sol);
53  v=1/5.0*RHS;
54  auto linferror = LInfError(v, sol);
55 
56  /* std::cout << "L2 error: " << l2error << std::endl;
57  std::cout << "L∞ Error: " << linferror << std::endl;*/
58  BOOST_REQUIRE(linferror < 1e-5);
59 
60  //domain.write("FDSOLVER_Lap_test");
61 
62  //domain.write("basic_test");
63 }
64 
65 
66  BOOST_AUTO_TEST_CASE(solver_Lap)
67  {
68  const size_t sz[2] = {82,82};
69  Box<2, double> box({0, 0}, {1, 1});
70  periodicity<2> bc = {NON_PERIODIC, NON_PERIODIC};
71  Ghost<2,long int> ghost(1);
72 
73  grid_dist_id<2, double, aggregate<double,double,double>> domain(sz, box, ghost, bc);
74 
75 
76  auto it = domain.getDomainIterator();
77  while (it.isNext())
78  {
79  auto key = it.get();
80  auto gkey = it.getGKey(key);
81  double x = gkey.get(0) * domain.spacing(0);
82  double y = gkey.get(1) * domain.spacing(1);
83  domain.get<0>(key) = sin(M_PI*x)*sin(M_PI*y);
84  domain.get<1>(key) = -2*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y);
85  ++it;
86  }
87 
88  domain.ghost_get<0>();
89  FD::Derivative_x Dx;
90  FD::Derivative_y Dy;
91  auto v = FD::getV<0>(domain);
92  auto RHS= FD::getV<1>(domain);
93  auto sol= FD::getV<2>(domain);
94 
95 
96  FD_scheme<equations2d1,decltype(domain)> Solver(ghost,domain);
97  FD::Lap Lap;
98  FD::LInfError LInfError;
99 
100 /* Solver.impose(Lap(v),{1,1},{79,79}, prop_id<1>());
101  Solver.impose(v,{0,0},{80,0}, prop_id<0>());
102  Solver.impose(v,{0,1},{0,79}, prop_id<0>());
103  Solver.impose(v,{0,80},{80,80}, prop_id<0>());
104  Solver.impose(v,{80,1},{80,79}, prop_id<0>());
105  Solver.solve(sol);*/
106 
107  Solver.impose(Lap(v),{1,1},{80,80}, prop_id<1>());
108  Solver.impose(v,{0,0},{81,0}, prop_id<0>());
109  Solver.impose(v,{0,1},{0,80}, prop_id<0>());
110  Solver.impose(v,{0,81},{81,81}, prop_id<0>());
111  Solver.impose(v,{81,1},{81,80}, prop_id<0>());
112  /*auto A=Solver.getA();
113  A.write("Lap_Matrix");*/
114  petsc_solver<double> pet_sol;
115  pet_sol.setPreconditioner(PCNONE);
116  Solver.solve_with_solver(pet_sol,sol);
117 
118  auto linferror = LInfError(v, sol);
119  BOOST_REQUIRE(linferror < 1e-3);
120  //std::cout << "L∞ Error: " << linferror << std::endl;
121 
122  //domain.write("FDSOLVER_Lap_test");
123  }
124 
125  BOOST_AUTO_TEST_CASE(solver_Lap_stag)
126  {
127  const size_t sz[2] = {82,82};
128  Box<2, double> box({0, 0}, {1, 1});
129  periodicity<2> bc = {NON_PERIODIC, NON_PERIODIC};
130  Ghost<2,long int> ghost(1);
131 
132  staggered_grid_dist<2, double, aggregate<double,double,double>> domain(sz, box, ghost, bc);
133 
134 
135  auto it = domain.getDomainIterator();
136  while (it.isNext())
137  {
138  auto key = it.get();
139  auto gkey = it.getGKey(key);
140  double x = gkey.get(0) * domain.spacing(0);
141  double y = gkey.get(1) * domain.spacing(1);
142  domain.get<0>(key) = sin(M_PI*x)*sin(M_PI*y);
143  domain.get<1>(key) = -2*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y);
144  ++it;
145  }
146 
147  domain.ghost_get<0>();
148  FD::Derivative_x Dx;
149  FD::Derivative_y Dy;
150  auto v = FD::getV<0>(domain);
151  auto RHS= FD::getV<1>(domain);
152  auto sol= FD::getV<2>(domain);
153 
154 
155  FD_scheme<equations2d1_stag,decltype(domain)> Solver(ghost,domain);
156  FD::Lap Lap;
157 
158 /* Solver.impose(Lap(v),{1,1},{79,79}, prop_id<1>());
159  Solver.impose(v,{0,0},{80,0}, prop_id<0>());
160  Solver.impose(v,{0,1},{0,79}, prop_id<0>());
161  Solver.impose(v,{0,80},{80,80}, prop_id<0>());
162  Solver.impose(v,{80,1},{80,79}, prop_id<0>());
163  Solver.solve(sol);*/
164 
165  Solver.impose(Lap(v),{1,1},{80,80}, prop_id<1>());
166  Solver.impose(v,{0,0},{81,0}, prop_id<0>());
167  Solver.impose(v,{0,1},{0,80}, prop_id<0>());
168  Solver.impose(v,{0,81},{81,81}, prop_id<0>());
169  Solver.impose(v,{81,1},{81,80}, prop_id<0>());
170  /*auto A=Solver.getA();
171  A.write("Lap_Matrix");*/
172 
173  Solver.solve(sol);
174 
175  auto it2 = domain.getDomainIterator();
176  double worst = 0.0;
177  while (it2.isNext()) {
178  auto p = it2.get();
179 
180  if (fabs(domain.getProp<0>(p) - domain.getProp<2>(p)) > worst) {
181  worst = fabs(domain.getProp<0>(p) - domain.getProp<2>(p));
182  }
183 
184  ++it2;
185  }
186  BOOST_REQUIRE(worst < 1e-3);
187  //std::cout << "Maximum Error: " << worst << std::endl;
188  //domain.write("FDSOLVER_Lap_test");
189  }
190  // Test failing for cores>=3
191  /*BOOST_AUTO_TEST_CASE(Lid_driven_PC)
192  { using namespace FD;
193  timer tt2;
194  tt2.start();
195  size_t gd_sz = 81;
196  double Re=100;
197  double V_err_eps = 1e-2;
198  double alpha=0.01;
199  constexpr int x = 0;
200  constexpr int y = 1;
201  const size_t szu[2] = {gd_sz,gd_sz};
202  int sz[2]={int(gd_sz),int(gd_sz)};
203  Box<2, double> box({0, 0}, {1,1});
204  periodicity<2> bc = {NON_PERIODIC, NON_PERIODIC};
205  double spacing;
206  spacing = 1.0 / (sz[0] - 1);
207  Ghost<2,long int> ghost(1);
208  auto &v_cl = create_vcluster();
209  //szu[0] = (size_t)sz[0];
210  //szu[1] = (size_t)sz[1];
211  typedef aggregate<double, VectorS<2, double>, VectorS<2, double>,VectorS<2, double>,double,VectorS<2, double>,double,double> LidCavity;
212  grid_dist_id<2, double, LidCavity> domain(szu, box,ghost,bc);
213  double x0, y0, x1, y1;
214  x0 = box.getLow(0);
215  y0 = box.getLow(1);
216  x1 = box.getHigh(0);
217  y1 = box.getHigh(1);
218  auto it = domain.getDomainIterator();
219  while (it.isNext())
220  {
221  auto key = it.get();
222  auto gkey = it.getGKey(key);
223  double x = gkey.get(0) * domain.spacing(0);
224  double y = gkey.get(1) * domain.spacing(1);
225  if (y==1)
226  {domain.get<1>(key) = 1.0;}
227  else
228  {domain.get<1>(key) = 0.0;}
229 
230  ++it;
231  }
232  domain.ghost_get<0>();
233  Derivative_x Dx;
234  Derivative_y Dy;
235  Derivative_xx Dxx;
236  Derivative_yy Dyy;
237  auto P = getV<0>(domain);
238  auto V = getV<1>(domain);
239  auto RHS = getV<2>(domain);
240  auto dV = getV<3>(domain);
241  auto div = getV<4>(domain);
242  auto V_star = getV<5>(domain);
243  auto H = getV<6>(domain);
244  auto dP = getV<7>(domain);
245 
246  //0 doesnt work
247  P = 0.0;
248  int n = 0, nmax = 50, ctr = 0, errctr=1, Vreset = 0;
249  double V_err=1;
250  if (Vreset == 1) {
251  Vreset = 0;
252  }
253  eq_id vx,vy;
254  vx.setId(0);
255  vy.setId(1);
256 
257  double sum, sum1, sum_k,V_err_old;
258  auto StokesX=(V[x]*Dx(V_star[x])+V[y]*Dy(V_star[x]))-(1.0/Re)*(Dxx(V_star[x])+Dyy(V_star[x]));
259  auto StokesY=(V[x]*Dx(V_star[y])+V[y]*Dy(V_star[y]))-(1.0/Re)*(Dxx(V_star[y])+Dyy(V_star[y]));
260  petsc_solver<double> solverPetsc;
261  solverPetsc.setSolver(KSPGMRES);
262 
263  RHS=0;
264  dV=0;
265 
266  timer tt;
267  while (V_err >= V_err_eps && n <= nmax) {
268  if (n%5==0){
269  domain.ghost_get<0,1 >(SKIP_LABELLING);
270  domain.write_frame("LID",n,BINARY);
271  domain.ghost_get<0>();
272  }
273  tt.start();
274  domain.ghost_get<0>(SKIP_LABELLING);
275  RHS[x] = -Dx(P);
276  RHS[y] = -Dy(P);
277  FD_scheme<equations2d2,decltype(domain)> Solver(ghost,domain);
278  Solver.impose(StokesX, {1,1},{sz[0]-2,sz[1]-2}, RHS[x], vx);
279  Solver.impose(StokesY, {1,1},{sz[0]-2,sz[1]-2}, RHS[y], vy);
280  Solver.impose(V_star[x], {0,sz[1]-1},{sz[0]-1,sz[1]-1}, 1.0, vx);
281  Solver.impose(V_star[y], {0,sz[1]-1},{sz[0]-1,sz[1]-1}, 0.0, vy);
282  Solver.impose(V_star[x], {0,1},{0,sz[1]-2}, 0.0, vx);
283  Solver.impose(V_star[y], {0,1},{0,sz[1]-2}, 0.0, vy);
284  Solver.impose(V_star[x], {sz[0]-1,1},{sz[0]-1,sz[1]-2}, 0.0, vx);
285  Solver.impose(V_star[y], {sz[0]-1,1},{sz[0]-1,sz[1]-2}, 0.0, vy);
286  Solver.impose(V_star[x], {0,0},{sz[0]-1,0}, 0.0, vx);
287  Solver.impose(V_star[y], {0,0},{sz[0]-1,0}, 0.0, vy);
288  //auto A=Solver.getA();
289  //A.write("Matrix");
290  Solver.solve(V_star[x], V_star[y]);
291  //return;
292  //Solver.solve_with_solver(solverPetsc,V_star[x], V_star[y]);
293 
294  domain.ghost_get<5>(SKIP_LABELLING);
295  div = (Dx(V_star[x]) + Dy(V_star[y]));
296 
297  FD_scheme<equations2d1E,decltype(domain)> SolverH(ghost,domain);
298  auto Helmholtz = Dxx(H)+Dyy(H);
299  SolverH.impose(Helmholtz,{1,1},{sz[0]-2,sz[1]-2},div);
300  SolverH.impose(Dy(H), {0,sz[0]-1},{sz[0]-1,sz[1]-1},0);
301  SolverH.impose(Dx(H), {sz[0]-1,1},{sz[0]-1,sz[1]-2},0);
302  SolverH.impose(Dx(H), {0,0},{sz[0]-1,0},0,vx,true);
303  SolverH.impose(H, {0,0},{0,0},0);
304  SolverH.impose(Dy(H), {0,1},{0,sz[1]-2},0);
305  //SolverH.solve_with_solver(solverPetsc2,H);
306  SolverH.solve(H);
307  //dP_bulk=Bulk_Lap(H);
308  domain.ghost_get<1,4,6>(SKIP_LABELLING);
309  P = P - alpha*(div-0.5*(V[x]*Dx(H)+V[y]*Dy(H)));
310  //dV[x]=Dx(H);
311  //dV[y]=Dy(H);
312  V_star[0] = V_star[0] - Dx(H);
313  V_star[1] = V_star[1] - Dy(H);
314  sum = 0;
315  sum1 = 0;
316  auto it2 = domain.getDomainIterator();
317  while (it2.isNext()) {
318  auto p = it2.get();
319  sum += (domain.getProp<5>(p)[0] - domain.getProp<1>(p)[0]) *
320  (domain.getProp<5>(p)[0] - domain.getProp<1>(p)[0]) +
321  (domain.getProp<5>(p)[1] - domain.getProp<1>(p)[1]) *
322  (domain.getProp<5>(p)[1] - domain.getProp<1>(p)[1]);
323  sum1 += domain.getProp<5>(p)[0] * domain.getProp<5>(p)[0] +
324  domain.getProp<5>(p)[1] * domain.getProp<5>(p)[1];
325  ++it2;
326  }
327 
328  V[x] = V_star[x];
329  V[y] = V_star[y];
330  v_cl.sum(sum);
331  v_cl.sum(sum1);
332  v_cl.execute();
333  sum = sqrt(sum);
334  sum1 = sqrt(sum1);
335  V_err_old = V_err;
336  V_err = sum / sum1;
337  if (V_err > V_err_old || abs(V_err_old - V_err) < 1e-8) {
338  errctr++;
339  //alpha_P -= 0.1;
340  } else {
341  errctr = 0;
342  }
343  if (n > 3) {
344  if (errctr > 5) {
345  //std::cout << "CONVERGENCE LOOP BROKEN DUE TO INCREASE/VERY SLOW DECREASE IN ERROR" << std::endl;
346  std::cout << "Alpha Halfed DUE TO INCREASE/VERY SLOW DECREASE IN ERROR" << std::endl;
347  //alpha=alpha/2;
348  Vreset = 1;
349  errctr=0;
350  //break;
351  } else {
352  Vreset = 0;
353  }
354  }
355  n++;
356  tt.stop();
357  if (v_cl.rank() == 0) {
358  std::cout << "Rel l2 cgs err in V = " << V_err << " at " << n << " and took " <<tt.getwct() <<"("<<tt.getcputime()<<") seconds(CPU)." << std::endl;
359  }
360  }
361  double worst1 = 0.0;
362  double worst2 = 0.0;
363 
364  domain.write("LID_final");
365  tt2.stop();
366  if (v_cl.rank() == 0) {
367  std::cout << "The simulation took " << tt2.getcputime() << "(CPU) ------ " << tt2.getwct()
368  << "(Wall) Seconds.";
369  }
370  }*/
371 
372 
373 
374 
375 
376 
377  /*
378 In 3D we use exact solution:
379 
380 Vx = x^2 + y^2
381 Vy = y^2 + z^2
382 Vz = x^2 + y^2 - 2(x+y)z
383 p = x + y + z - 3/2
384 f_x = f_y = f_z = 3
385 
386 -\Delta u + \nabla p + f = <-4, -4, -4> + <1, 1, 1> + <3, 3, 3> = 0
387 \nabla \cdot u = 2x + 2y - 2(x + y) = 0
388 */
389 
390  struct lid_nn_stag
391  {
392  // dimensionaly of the equation (2D problem 3D problem ...)
393  static const unsigned int dims = 2;
394 
395  // number of fields in the system v_x, v_y, P so a total of 3
396  static const unsigned int nvar = 3;
397 
398  // boundary conditions PERIODIC OR NON_PERIODIC
399  static const bool boundary[];
400 
401  // type of space float, double, ...
402  typedef float stype;
403 
404  // type of base grid, it is the distributed grid that will store the result
405  // Note the first property is a 2D vector (velocity), the second is a scalar (Pressure)
406  typedef staggered_grid_dist<2,float,aggregate<float[2],float>,CartDecomposition<2,float>> b_grid;
407 
408  // type of SparseMatrix, for the linear system, this parameter is bounded by the solver
409  // that you are using, in case of umfpack it is the only possible choice
410  typedef SparseMatrix<double,int> SparseMatrix_type;
411 
412  // type of Vector for the linear system, this parameter is bounded by the solver
413  // that you are using, in case of umfpack it is the only possible choice
414  typedef Vector<double> Vector_type;
415 
416  typedef umfpack_solver<double> solver_type;
417  };
418 
419  const bool lid_nn_stag::boundary[] = {NON_PERIODIC,NON_PERIODIC};
420 
421  BOOST_AUTO_TEST_CASE(solver_Lid_driven_cavity_stag)
422  {
423  Vcluster<> & v_cl = create_vcluster();
424 
425  if (v_cl.getProcessingUnits() > 3)
426  return;
427 
429 
430  // velocity in the grid is the property 0, pressure is the property 1
431  constexpr int velocity = 0;
432  constexpr int pressure = 1;
433 
434  constexpr int x = 0;
435  constexpr int y = 1;
436 
437  // Domain, a rectangle
438  Box<2,float> domain({0.0,0.0},{3.0,1.0});
439 
440  // Ghost (Not important in this case but required)
441  //Ghost<2,float> g(0.01); // suited for particles
442  Ghost<2,long int> g(1); //better for grids
443  // Grid points on x=256 and y=64
444  long int sz[] = {256,64};
445  size_t szu[2];
446  szu[0] = (size_t)sz[0];
447  szu[1] = (size_t)sz[1];
448 
449  // We need one more point on the left and down part of the domain
450  // This is given by the boundary conditions that we impose, the
451  // reason is mathematical in order to have a well defined system
452  // and cannot be discussed here
453  Padding<2> pd({1,1},{0,0});
454 
455  // Distributed grid that store the solution
456  staggered_grid_dist<2,float,aggregate<Point<2,float>,float,Point<2,float>,float>> g_dist(szu,domain,g);
457  grid_dist_id<2,float,aggregate<Point<2,float>,float,Point<2,float>,float>> g_dist_normal(g_dist.getDecomposition(),szu,g);
458 
459  openfpm::vector<comb<2>> cmb_v;
460  cmb_v.add({0,-1});
461  cmb_v.add({-1,0});
462 
463  g_dist.setDefaultStagPosition();
464  g_dist.setStagPosition<0>(cmb_v);
465 
466  // It is the maximum extension of the stencil
467  Ghost<2,long int> stencil_max(1);
468 
469  // Finite difference scheme
470  FD_scheme<lid_nn_stag,decltype(g_dist)> fd(pd,stencil_max,g_dist);
471 
472  auto P = FD::getV_stag<1>(g_dist);
473  auto v = FD::getV_stag<0>(g_dist);
474  auto RHS = FD::getV_stag<2>(g_dist);
475  //auto RHSy = FD::getV_stag<3>(g_dist);
476 
477 
478  v.setVarId(0);
479  P.setVarId(2);
480 
481  FD::Derivative_x_stag Dx;
482  FD::Derivative_y_stag Dy;
483  FD::Lap Lap;
484 
485  double nu = 1.0;
486 
487  auto Stokes_vx = nu * Lap(v[x]) - Dx(P);
488  auto Stokes_vy = nu * Lap(v[y]) - Dy(P);
489 
490  auto incompressibility = Dx(v[x]) + Dy(v[y]);
491 
492  eq_id ic,vx,vy;
493 
494  ic.setId(2);
495  vx.setId(0);
496  vy.setId(1);
497 
498  comb<2> center_cell({0,0});
499  comb<2> left_cell({0,-1});
500  comb<2> bottom_cell({-1,0});
501  comb<2> corner_right({-1,1});
502  comb<2> corner_dw({-1,-1});
503  comb<2> corner_up({1,-1});
504 
505  auto it = g_dist.getDomainIterator();
506  while (it.isNext())
507  {
508  auto key = it.get();
509  auto gkey = it.getGKey(key);
510  double xg = gkey.get(0) * g_dist.spacing(0);
511  double yg = gkey.get(1) * g_dist.spacing(1);
512  if (xg==3.0){
513  //v[y].value(key,bottom_cell)=1.0;
514  g_dist.getProp<2>(key)[1] = 1.0;
515  g_dist.getProp<3>(key) = 1.0;
516  //std::cout<<it.getGKey(key).get(0)<<","<<it.getGKey(key).get(1)<<":"<<g_dist.getProp<3>(key)<<std::endl;
517  }
518  ++it;
519  }
520 
521  // Here we impose the equation, we start from the incompressibility Eq imposed in the bulk with the
522  // exception of the first point {0,0} and than we set P = 0 in {0,0}, why we are doing this is again
523  // mathematical to have a well defined system, an intuitive explanation is that P and P + c are both
524  // solution for the incompressibility equation, this produce an ill-posed problem to make it well posed
525  // we set one point in this case {0,0} the pressure to a fixed constant for convenience P = 0
526 
527  fd.impose(incompressibility, {0,0},{sz[0]-2,sz[1]-2}, 0.0,ic,true);
528  fd.impose(P, {0,0},{0,0},0.0,ic);
529 
530  // Here we impose the Eq1 and Eq2
531  fd.impose(Stokes_vx, {1,0},{sz[0]-2,sz[1]-2},0.0,vx,left_cell);
532  fd.impose(Stokes_vy, {0,1},{sz[0]-2,sz[1]-2},0.0,vy,bottom_cell);
533 
534  // Staggering pattern in the domain
535  // +--Vy-+
536  // | |
537  // Vx P Vx
538  // | |
539  // 0--Vy-+
540  //
541  // Imposing B1
542  // Left Wall Vx without touching top wall and with specified location in the cell ("left_cell").
543  fd.impose(v[x], {0,0},{0,sz[1]-2},0.0,vx,left_cell);
544  // Left Wall Vy without touching top wall
545  // (Starts at -1 index since Vy is on top of the cell and needs special boundary treatment "corner_dw")
546  // cells in the Domain near 0,0 look like:
547  //
548  // | |
549  // Vx P(0,0)
550  // | |
551  // Vy---Vy--+
552  // :
553  // "corner_right of -1" specified by starting at -1,0
554  fd.impose(v[y], {-1,0},{-1,sz[1]-1},0.0,vy,corner_right);
555 
556  //Imposing B2
557  // Similarly Right Wall (Also need "corner_dw" treatment for Vy, hence +1 in the y index)
558  fd.impose(v[x],{sz[0]-1,0},{sz[0]-1,sz[1]-2},0.0,vx,left_cell);
559  fd.impose(v[y],{sz[0]-1,0},{sz[0]-1,sz[1]-1},prop_id<3>(),vy,corner_dw);
560 
561  // Imposing B3
562  // Similarly Bottom Wall (needs "corner_up" treatment for Vx, hence -1 in the x index)
563  fd.impose(v[x], {0,-1},{sz[0]-1,-1},0.0,vx,corner_up);
564  fd.impose(v[y], {0,0},{sz[0]-2,0},0.0,vy,bottom_cell);
565  // Imposing B4
566  // Similarly Top Wall (needs "corner_dw" treatment for Vx, hence +1 in the x index)
567  fd.impose(v[x],{0,sz[1]-1},{sz[0]-1,sz[1]-1},0.0,vx,corner_dw);
568  fd.impose(v[y],{0,sz[1]-1},{sz[0]-2,sz[1]-1},0.0,vy,bottom_cell);
569 
570  // When we pad the grid, there are points of the grid that are not
571  // touched by the previous condition. Mathematically this lead
572  // to have too many variables for the conditions that we are imposing.
573  // Here we are imposing variables that we do not touch to zero
574  //
575 
576  // Padding pressure
577  fd.impose(P, {-1,-1},{sz[0]-1,-1},0.0,ic);
578  fd.impose(P, {-1,sz[1]-1},{sz[0]-1,sz[1]-1},0.0,ic);
579  fd.impose(P, {-1,0},{-1,sz[1]-2},0.0,ic);
580  fd.impose(P, {sz[0]-1,0},{sz[0]-1,sz[1]-2},0.0,ic);
581 
582  // Impose v_x Padding Impose v_y padding
583  fd.impose(v[x], {-1,-1},{-1,sz[1]-1},0.0,vx,left_cell);
584  fd.impose(v[y], {-1,-1},{sz[0]-1,-1},0.0,vy,bottom_cell);
585 
586  fd.solve(v[x],v[y],P);
587 
588  // auto it5 = g_dist.getSubDomainIterator({1,1},{(int)g_dist.size(0)-1,(int)g_dist.size(1)-1});
589 
590  // while (it5.isNext())
591  // {
592  // auto key2 = it5.get();
593 
594  // std::cout << g_dist.template getProp<0>(key2)[0] << " " << std::endl;
595 
596  // ++it5;
597  // }
598 
599  g_dist.template ghost_get<0,1>();
600 
601  // Carefull we cannot interpolate on the border
602 
603  auto it3 = g_dist_normal.getSubDomainIterator({1,1},{(int)g_dist_normal.size(0)-1,(int)g_dist_normal.size(1)-1});
604  auto it4 = g_dist.getSubDomainIterator({1,1},{(int)g_dist.size(0)-1,(int)g_dist.size(1)-1});
605 
606  while (it3.isNext())
607  {
608  auto key = it3.get();
609  auto key2 = it4.get();
610 
611  g_dist_normal.template getProp<0>(key)[0] = v[x].value(key2,corner_dw);
612  g_dist_normal.template getProp<0>(key)[1] = v[y].value(key2,corner_dw);;
613 
614  g_dist_normal.template getProp<1>(key) = P.value(key2,corner_dw);
615  g_dist_normal.template getProp<2>(key)[x] = RHS[x].value(key2,corner_dw);
616  g_dist_normal.template getProp<2>(key)[y] = RHS[y].value(key2,corner_dw);
617 
618 
619  ++it3;
620  ++it4;
621  }
622 
623 #if !(defined(SE_CLASS3) || defined(TEST_COVERAGE_MODE))
624 
625  // Initialize openfpm
626  grid_dist_id<2,float,aggregate<float[2],float>> g_dist2(g_dist.getDecomposition(),szu,g);
627  g_dist2.load("test/lid_driven_cavity_reference.hdf5");
628 
629  auto it2 = g_dist2.getSubDomainIterator({1,1},{(int)g_dist2.size(0)-1,(int)g_dist2.size(1)-1});
630 
631  bool test = true;
632  while (it2.isNext())
633  {
634  auto p = it2.get();
635 
636  test &= fabs(g_dist2.template getProp<velocity>(p)[0] - g_dist_normal.template getProp<velocity>(p)[0]) < 3.5e-3;
637  test &= fabs(g_dist2.template getProp<velocity>(p)[1] - g_dist_normal.template getProp<velocity>(p)[1]) < 3.5e-3;
638 
639  test &= fabs(g_dist2.template getProp<pressure>(p) - g_dist_normal.template getProp<pressure>(p)) < 3.0;
640 
641  if (test == false)
642  {
643  std::cout << g_dist2.template getProp<velocity>(p)[0] << " " << g_dist_normal.template getProp<velocity>(p)[0] << std::endl;
644  std::cout << g_dist2.template getProp<velocity>(p)[1] << " " << g_dist_normal.template getProp<velocity>(p)[1] << std::endl;
645 
646  std::cout << g_dist2.template getProp<pressure>(p) << " " << g_dist_normal.template getProp<pressure>(p) << std::endl;
647 
648  //g_dist_normal.template getProp<0>(p)[1] = v[y].value(p,corner_dw);
649 
650  g_dist_normal.template getProp<1>(p) = P.value(p,corner_dw);
651 
652  break;
653  }
654 
655  ++it2;
656  }
657 
658  BOOST_REQUIRE_EQUAL(test,true);
659 
660 #endif
661  }
662 
663 
664  struct ana_nn_stag
665  {
666  // dimensionaly of the equation (2D problem 3D problem ...)
667  static const unsigned int dims = 2;
668 
669  // number of fields in the system v_x, v_y, P so a total of 3
670  static const unsigned int nvar = 3;
671 
672  // boundary conditions PERIODIC OR NON_PERIODIC
673  static const bool boundary[];
674 
675  // type of space float, double, ...
676  typedef float stype;
677 
678  // type of base grid, it is the distributed grid that will store the result
679  // Note the first property is a 2D vector (velocity), the second is a scalar (Pressure)
680  typedef staggered_grid_dist<2,float,aggregate<float[2],float,float[2],float>,CartDecomposition<2,float>> b_grid;
681 
682  // type of SparseMatrix, for the linear system, this parameter is bounded by the solver
683  // that you are using, in case of umfpack it is the only possible choice
684  // typedef SparseMatrix<double,int> SparseMatrix_type;
685  // // typedef SparseMatrix<double, int, PETSC_BASE> SparseMatrix_type;
686 
687  // // type of Vector for the linear system, this parameter is bounded by the solver
688  // // that you are using, in case of umfpack it is the only possible choice
689  // typedef Vector<double> Vector_type;
690  // // typedef Vector<double, PETSC_BASE> Vector_type;
691 
692  // typedef umfpack_solver<double> solver_type;
693  // // typedef petsc_solver<double> solver_type;
695  typedef SparseMatrix<double, int, PETSC_BASE> SparseMatrix_type;
696 
698  typedef Vector<double, PETSC_BASE> Vector_type;
699 
700  typedef petsc_solver<double> solver_type;
701  };
702 
703  const bool ana_nn_stag::boundary[] = {NON_PERIODIC,NON_PERIODIC};
704 
705  /*
706  In 2D we use exact solution:
707 
708  u = x^2 + y^2
709  v = 2 x^2 - 2xy
710  p = x + y - 1
711  f_x = f_y = 3
712 
713  so that
714 
715  -\Delta u + \nabla p + f = <-4, -4> + <1, 1> + <3, 3> = 0
716  \Delta u - \nabla p = +f
717  \nabla \cdot u = 2x - 2x = 0
718  */
719 #if 0
720  //These tests need development and checks
721  BOOST_AUTO_TEST_CASE(solver_ana_stag)
722  {
723  Vcluster<> & v_cl = create_vcluster();
724 
725  if (v_cl.getProcessingUnits() > 3)
726  return;
727 
728  // velocity in the grid is the property 0, pressure is the property 1
729  constexpr int velocity = 0;
730  constexpr int pressure = 1;
731 
732  constexpr int x = 0;
733  constexpr int y = 1;
734 
735  // Domain, a rectangle
736  Box<2,float> domain({0.0,0.0},{1.0,1.0});
737 
738  // Grid points on x=256 and y=256
739  long int sz[] = {256,256};
740  // Ghost (Not important in this case but required)
741  Ghost<2,float> g(domain.getHigh(0)/sz[0]);
742 
743  size_t szu[2];
744  szu[0] = (size_t)sz[0];
745  szu[1] = (size_t)sz[1];
746 
747  Padding<2> pd({1,1},{0,0});
748 
749  // Distributed grid that store the solution
750  staggered_grid_dist<2,float,aggregate<Point<2,float>,float,Point<2,float>,float,float,float>> g_dist(szu,domain,g);
751  grid_dist_id<2,float,aggregate<Point<2,float>,float,Point<2,float>,float,float,float>> g_dist_normal(g_dist.getDecomposition(),szu,g);
752 
753  double hx = g_dist.spacing(0);
754  double hy = g_dist.spacing(1);
755 
756  openfpm::vector<comb<2>> cmb_v;
757  cmb_v.add({0,-1});
758  cmb_v.add({-1,0});
759 
760  g_dist.setDefaultStagPosition();
761  g_dist.setStagPosition<0>(cmb_v);
762 
763  // It is the maximum extension of the stencil
764  Ghost<2,long int> stencil_max(1);
765 
766  // Finite difference scheme
767  FD_scheme<lid_nn_stag,decltype(g_dist)> fd(pd,stencil_max,g_dist);
768 
769  auto P = FD::getV_stag<1>(g_dist);
770  auto v = FD::getV_stag<0>(g_dist);
771  auto Ana_v = FD::getV_stag<2>(g_dist);
772  auto Ana_P = FD::getV_stag<3>(g_dist);
773 
774  v.setVarId(0);
775  P.setVarId(2);
776 
777  FD::Derivative_x_stag Dx;
778  FD::Derivative_y_stag Dy;
779  FD::Lap Lap;
780 
781  auto Stokes_vx = Lap(v[x]) - Dx(P);
782  auto Stokes_vy = Lap(v[y]) - Dy(P);
783 
784  auto incompressibility = Dx(v[x]) + Dy(v[y]);
785 
786  eq_id ic,vx,vy;
787 
788  ic.setId(2);
789  vx.setId(0);
790  vy.setId(1);
791 
792  comb<2> center_cell({0,0});
793  comb<2> left_cell({0,-1});
794  comb<2> bottom_cell({-1,0});
795  comb<2> corner_right({-1,1});
796  comb<2> corner_dw({-1,-1});
797  comb<2> corner_up({1,-1});
798 
799  auto fu = [](double x, double y) { return x*x + y*y; };
800  auto fv = [](double x, double y) { return 2.0*x*x - 2.0*x*y; };
801 
802  // Staggering pattern in the domain
803  // +--Vy-+
804  // | |
805  // Vx P Vx
806  // | |
807  // 0--Vy-+
808  //
809 
810  fd.impose(incompressibility, {0,0},{sz[0]-2,sz[1]-2}, 0.0,ic,true);
811  fd.impose(P, {0,0},{0,0},(hx+hy)/2-1.0,ic);
812 
813  fd.impose(Stokes_vx, {1,0},{sz[0]-2,sz[1]-2},3.0,vx,left_cell);
814  fd.impose(Stokes_vy, {0,1},{sz[0]-2,sz[1]-2},3.0,vy,bottom_cell);
815 
816  // Imposing B1 Left Wall
817  fd.impose(v[x], {0,0},{0,sz[1]-2},fu,vx,left_cell);
818  fd.impose(v[y], {-1,0},{-1,sz[1]-1},fv,vy,corner_right);
819 
820  // //Imposing B2 Right Wall
821  fd.impose(v[x],{sz[0]-1,0},{sz[0]-1,sz[1]-2},fu,vx,left_cell);
822  fd.impose(v[y],{sz[0]-1,0},{sz[0]-1,sz[1]-1},fv,vy,corner_dw);
823 
824  // Imposing B3 Bottom Wall
825  fd.impose(v[x],{0,-1},{sz[0]-1,-1},fu,vx,corner_up);
826  fd.impose(v[y], {0,0},{sz[0]-2,0},fv,vy,bottom_cell);
827 
828  // Imposing B4 Top Wall
829  fd.impose(v[x],{0,sz[1]-1},{sz[0]-1,sz[1]-1},fu,vx,corner_dw);
830  fd.impose(v[y],{0,sz[1]-1},{sz[0]-2,sz[1]-1},fv,vy,bottom_cell);
831 
832  // Padding pressure
833  fd.impose(P, {-1,-1},{sz[0]-1,-1},0.0,ic);
834  fd.impose(P, {-1,sz[1]-1},{sz[0]-1,sz[1]-1},0.0,ic);
835  fd.impose(P, {-1,0},{-1,sz[1]-2},0.0,ic);
836  fd.impose(P, {sz[0]-1,0},{sz[0]-1,sz[1]-2},0.0,ic);
837 
838  // Impose v_x Padding Impose v_y padding
839  fd.impose(v[x], {-1,-1},{-1,sz[1]-1},0.0,vx,left_cell);
840  fd.impose(v[y], {-1,-1},{sz[0]-1,-1},0.0,vy,bottom_cell);
841 
842  fd.solve(v[x],v[y],P);
843 
844  auto it2 = g_dist_normal.getDomainIterator();
845 
846  while (it2.isNext())
847  {
848  auto key = it2.get();
849  auto gkey = it2.getGKey(key);
850  double x = gkey.get(0) * g_dist_normal.spacing(0);
851  double y = gkey.get(1) * g_dist_normal.spacing(1);
852 
853  g_dist_normal.template getProp<0>(key)[0] = v[0].value(key,left_cell);
854  g_dist_normal.template getProp<0>(key)[1] = v[1].value(key,bottom_cell);
855 
856  g_dist_normal.template getProp<1>(key) = P.value(key,center_cell);
857 
858  g_dist_normal.template getProp<2>(key)[0] = x*x + y*y;
859  g_dist_normal.template getProp<2>(key)[1] = 2*x*x - 2*x*y;
860  g_dist_normal.template getProp<3>(key) = x + y - 1;
861  ++it2;
862  }
863 
864  g_dist_normal.write("ana_stokes");
865 
867 
868  //g_dist.write("ana_stokes")
869  }
870 
871 #endif
872 
873 
874 BOOST_AUTO_TEST_SUITE_END()
875 #endif //Have Eigen
876 #endif //Have Petsc
Box
This class represent an N-dimensional box.
Definition: Box.hpp:60
CartDecomposition
This class decompose a space into sub-sub-domains and distribute them across processors.
Definition: CartDecomposition.hpp:145
FD_scheme
Finite Differences.
Definition: FD_Solver.hpp:128
Ghost
Definition: Ghost.hpp:40
Lap
Laplacian second order on h (spacing)
Definition: Laplacian.hpp:23
Padding
Class that contain Padding information on each direction positive and Negative direction.
Definition: Ghost.hpp:127
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
SparseMatrix
Sparse Matrix implementation.
Definition: SparseMatrix.hpp:79
Vcluster_base::getProcessingUnits
size_t getProcessingUnits()
Get the total number of processors.
Definition: VCluster_base.hpp:501
Vcluster
Implementation of VCluster class.
Definition: VCluster.hpp:59
Vector
Sparse Matrix implementation stub object when OpenFPM is compiled with no linear algebra support.
Definition: Vector.hpp:40
eq_id
Definition: eq_solve_common.hpp:15
grid_dist_id
This is a distributed grid.
Definition: grid_dist_id.hpp:279
openfpm::vector
Implementation of 1-D std::vector like structure.
Definition: map_vector.hpp:204
petsc_solver
This class is able to do Matrix inversion in parallel with PETSC solvers.
Definition: petsc_solver.hpp:102
petsc_solver::setPreconditioner
void setPreconditioner(PCType type)
Set the preconditioner of the linear solver.
Definition: petsc_solver.hpp:1046
staggered_grid_dist
Implementation of the staggered grid.
Definition: staggered_dist_grid.hpp:64
umfpack_solver< double >
stub when library compiled without eigen
Definition: umfpack_solver.hpp:219
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
comb
Position of the element of dimension d in the hyper-cube of dimension dim.
Definition: comb.hpp:35
equations2d1_stag
Specify the general characteristic of system to solve.
Definition: EqnsStructFD.hpp:478
equations2d1
Specify the general characteristic of system to solve.
Definition: EqnsStruct.hpp:14
periodicity
Boundary conditions.
Definition: common.hpp:22
prop_id
Definition: eq_solve_common.hpp:12