OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
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
17BOOST_AUTO_TEST_SUITE( FD_Solver_test )
18
19
20BOOST_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 /*
378In 3D we use exact solution:
379
380Vx = x^2 + y^2
381Vy = y^2 + z^2
382Vz = x^2 + y^2 - 2(x+y)z
383p = x + y + z - 3/2
384f_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
874BOOST_AUTO_TEST_SUITE_END()
875#endif //Have Eigen
876#endif //Have Petsc
Box
This class represent an N-dimensional box.
Definition Box.hpp:61
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
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:479
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:278
openfpm::vector
Implementation of 1-D std::vector like structure.
Definition map_vector.hpp:203
petsc_solver
In case T does not match the PETSC precision compilation create a stub structure.
Definition petsc_solver.hpp:56
staggered_grid_dist
Implementation of the staggered grid.
Definition staggered_dist_grid.hpp:64
umfpack_solver
stub when library compiled without eigen
Definition umfpack_solver.hpp:194
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