OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
main_petsc.cpp
1
33
34#include "config.h"
35
36#ifdef HAVE_PETSC
37
38#include "Grid/grid_dist_id.hpp"
39#include "Matrix/SparseMatrix.hpp"
40#include "Vector/Vector.hpp"
41#include "FiniteDifference/FDScheme.hpp"
42#include "FiniteDifference/util/common.hpp"
43#include "FiniteDifference/eq.hpp"
44#include "Solvers/umfpack_solver.hpp"
45#include "Solvers/petsc_solver.hpp"
46
47struct lid_nn
48{
49 // dimensionaly of the equation ( 3D problem ...)
50 static const unsigned int dims = 3;
51
52 // number of fields in the system v_x, v_y, v_z, P so a total of 4
53 static const unsigned int nvar = 4;
54
55 // boundary conditions PERIODIC OR NON_PERIODIC
56 static const bool boundary[];
57
58 // type of space float, double, ...
59 typedef float stype;
60
61 // type of base grid, it is the distributed grid that will store the result
62 // Note the first property is a 3D vector (velocity), the second is a scalar (Pressure)
64
65 // type of SparseMatrix, for the linear system, this parameter is bounded by the solver
66 // that you are using, in case of umfpack using <double,int> it is the only possible choice
67 // We select PETSC SparseMatrix implementation
68 typedef SparseMatrix<double,int,PETSC_BASE> SparseMatrix_type;
69
70 // type of Vector for the linear system, this parameter is bounded by the solver
71 // that you are using, in case of umfpack using <double> it is the only possible choice
72 // We select PETSC implementation
73 typedef Vector<double,PETSC_BASE> Vector_type;
74
75 // Define that the underline grid where we discretize the system of equation is staggered
76 static const int grid_type = STAGGERED_GRID;
77};
78
79const bool lid_nn::boundary[] = {NON_PERIODIC,NON_PERIODIC,NON_PERIODIC};
80
82
83
117
118// Constant Field
119struct eta
120{
122 typedef void const_field;
123
125 static float val() {return 1.0;}
126};
127
128// Model the equations
129
130constexpr unsigned int v[] = {0,1,2};
131constexpr unsigned int P = 3;
132constexpr unsigned int ic = 3;
133constexpr int x = 0;
134constexpr int y = 1;
135constexpr int z = 2;
136
137typedef Field<v[x],lid_nn> v_x;
138typedef Field<v[y],lid_nn> v_y;
139typedef Field<v[z],lid_nn> v_z;
140typedef Field<P,lid_nn> Prs;
141
142// Eq1 V_x
143
145typedef D<x,Prs,lid_nn> p_x; // Step 2
146typedef minus<p_x,lid_nn> m_p_x; // Step3
147typedef sum<eta_lap_vx,m_p_x,lid_nn> vx_eq; // Step4
148
149// Eq2 V_y
150
152typedef D<y,Prs,lid_nn> p_y;
155
156// Eq3 V_z
157
159typedef D<z,Prs,lid_nn> p_z;
162
163// Eq4 Incompressibility
164
165typedef D<x,v_x,lid_nn,FORWARD> dx_vx; // Step5
166typedef D<y,v_y,lid_nn,FORWARD> dy_vy; // Step6
167typedef D<z,v_z,lid_nn,FORWARD> dz_vz; // Step 7
168typedef sum<dx_vx,dy_vy,dz_vz,lid_nn> ic_eq; // Step8
169
171
206
207// Directional Avg
210
213
216
219
222
225
226
227#define EQ_1 0
228#define EQ_2 1
229#define EQ_3 2
230#define EQ_4 3
231
233
234#include "Vector/vector_dist.hpp"
235#include "data_type/aggregate.hpp"
236
237int main(int argc, char* argv[])
238{
259
260 // Initialize
261 openfpm_init(&argc,&argv);
262 {
263 // velocity in the grid is the property 0, pressure is the property 1
264 constexpr int velocity = 0;
265 constexpr int pressure = 1;
266
267 // Domain
268 Box<3,float> domain({0.0,0.0,0.0},{3.0,1.0,1.0});
269
270 // Ghost (Not important in this case but required)
271 Ghost<3,float> g(0.01);
272
273 // Grid points on x=36 and y=12 z=12
274 long int sz[] = {36,12,12};
275 size_t szu[3];
276 szu[0] = (size_t)sz[0];
277 szu[1] = (size_t)sz[1];
278 szu[2] = (size_t)sz[2];
279
280
281 // We need one more point on the left and down part of the domain
282 // This is given by the boundary conditions that we impose.
283 //
284 Padding<3> pd({1,1,1},{0,0,0});
285
287
300
302
304
323
324 // It is the maximum extension of the stencil (order 2 laplacian stencil has extension 1)
325 Ghost<3,long int> stencil_max(1);
326
327 // Finite difference scheme
328 FDScheme<lid_nn> fd(pd, stencil_max, domain, g_dist);
329
331
356
357 // start and end of the bulk
358
359 fd.impose(ic_eq(),0.0, EQ_4, {0,0,0},{sz[0]-2,sz[1]-2,sz[2]-2},true);
360 fd.impose(Prs(), 0.0, EQ_4, {0,0,0},{0,0,0});
361 fd.impose(vx_eq(),0.0, EQ_1, {1,0},{sz[0]-2,sz[1]-2,sz[2]-2});
362 fd.impose(vy_eq(),0.0, EQ_2, {0,1},{sz[0]-2,sz[1]-2,sz[2]-2});
363 fd.impose(vz_eq(),0.0, EQ_3, {0,0,1},{sz[0]-2,sz[1]-2,sz[2]-2});
364
365 // v_x
366 // R L (Right,Left)
367 fd.impose(v_x(),0.0, EQ_1, {0,0,0}, {0,sz[1]-2,sz[2]-2});
368 fd.impose(v_x(),0.0, EQ_1, {sz[0]-1,0,0},{sz[0]-1,sz[1]-2,sz[2]-2});
369
370 // T B (Top,Bottom)
371 fd.impose(avg_y_vx_f(),0.0, EQ_1, {0,-1,0}, {sz[0]-1,-1,sz[2]-2});
372 fd.impose(avg_y_vx(),0.0, EQ_1, {0,sz[1]-1,0},{sz[0]-1,sz[1]-1,sz[2]-2});
373
374 // A F (Forward,Backward)
375 fd.impose(avg_z_vx_f(),0.0, EQ_1, {0,-1,-1}, {sz[0]-1,sz[1]-1,-1});
376 fd.impose(avg_z_vx(),0.0, EQ_1, {0,-1,sz[2]-1},{sz[0]-1,sz[1]-1,sz[2]-1});
377
378 // v_y
379 // R L
380 fd.impose(avg_x_vy_f(),0.0, EQ_2, {-1,0,0}, {-1,sz[1]-1,sz[2]-2});
381 fd.impose(avg_x_vy(),1.0, EQ_2, {sz[0]-1,0,0},{sz[0]-1,sz[1]-1,sz[2]-2});
382
383 // T B
384 fd.impose(v_y(), 0.0, EQ_2, {0,0,0}, {sz[0]-2,0,sz[2]-2});
385 fd.impose(v_y(), 0.0, EQ_2, {0,sz[1]-1,0},{sz[0]-2,sz[1]-1,sz[2]-2});
386
387 // F A
388 fd.impose(avg_z_vy(),0.0, EQ_2, {-1,0,sz[2]-1}, {sz[0]-1,sz[1]-1,sz[2]-1});
389 fd.impose(avg_z_vy_f(),0.0, EQ_2, {-1,0,-1}, {sz[0]-1,sz[1]-1,-1});
390
391 // v_z
392 // R L
393 fd.impose(avg_x_vz_f(),0.0, EQ_3, {-1,0,0}, {-1,sz[1]-2,sz[2]-1});
394 fd.impose(avg_x_vz(),1.0, EQ_3, {sz[0]-1,0,0},{sz[0]-1,sz[1]-2,sz[2]-1});
395
396 // T B
397 fd.impose(avg_y_vz(),0.0, EQ_3, {-1,sz[1]-1,0},{sz[0]-1,sz[1]-1,sz[2]-1});
398 fd.impose(avg_y_vz_f(),0.0, EQ_3, {-1,-1,0}, {sz[0]-1,-1,sz[2]-1});
399
400 // F A
401 fd.impose(v_z(),0.0, EQ_3, {0,0,0}, {sz[0]-2,sz[1]-2,0});
402 fd.impose(v_z(),0.0, EQ_3, {0,0,sz[2]-1},{sz[0]-2,sz[1]-2,sz[2]-1});
403
404 // When we pad the grid, there are points of the grid that are not
405 // touched by the previous condition. Mathematically this lead
406 // to have too many variables for the conditions that we are imposing.
407 // Here we are imposing variables that we do not touch to zero
408 //
409
410 // L R
411 fd.impose(Prs(), 0.0, EQ_4, {-1,-1,-1},{-1,sz[1]-1,sz[2]-1});
412 fd.impose(Prs(), 0.0, EQ_4, {sz[0]-1,-1,-1},{sz[0]-1,sz[1]-1,sz[2]-1});
413
414 // T B
415 fd.impose(Prs(), 0.0, EQ_4, {0,sz[1]-1,-1}, {sz[0]-2,sz[1]-1,sz[2]-1});
416 fd.impose(Prs(), 0.0, EQ_4, {0,-1 ,-1}, {sz[0]-2,-1, sz[2]-1});
417
418 // F A
419 fd.impose(Prs(), 0.0, EQ_4, {0,0,sz[2]-1}, {sz[0]-2,sz[1]-2,sz[2]-1});
420 fd.impose(Prs(), 0.0, EQ_4, {0,0,-1}, {sz[0]-2,sz[1]-2,-1});
421
422 // Impose v_x v_y v_z padding
423 fd.impose(v_x(), 0.0, EQ_1, {-1,-1,-1},{-1,sz[1]-1,sz[2]-1});
424 fd.impose(v_y(), 0.0, EQ_2, {-1,-1,-1},{sz[0]-1,-1,sz[2]-1});
425 fd.impose(v_z(), 0.0, EQ_3, {-1,-1,-1},{sz[0]-1,sz[1]-1,-1});
426
428
448
449 // Create an PETSC solver
451
452 // Warning try many solver and collect statistics require a lot of time
453 // To just solve you can comment this line
454// solver.best_solve();
455
456 // Give to the solver A and b, return x, the solution
457 auto x = solver.solve(fd.getA(),fd.getB());
458
460
473
474 // copy the solution to grid
475 fd.template copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1,sz[2]-1},g_dist);
476
477 g_dist.write("lid_driven_cavity_p_petsc");
478
480
481
494 }
495 openfpm_finalize();
496
498
499
508}
509
510#else
511
512int main(int argc, char* argv[])
513{
514 return 0;
515}
516
517#endif
Average.
Definition Average.hpp:27
This class represent an N-dimensional box.
Definition Box.hpp:61
This class decompose a space into sub-sub-domains and distribute them across processors.
Derivative second order on h (spacing)
Finite Differences.
Definition FDScheme.hpp:127
Definition eq.hpp:83
Class that contain Padding information on each direction positive and Negative direction.
Definition Ghost.hpp:127
Test structure used for several test.
Sparse Matrix implementation.
Sparse Matrix implementation stub object when OpenFPM is compiled with no linear algebra support.
Definition Vector.hpp:40
This is a distributed grid.
In case T does not match the PETSC precision compilation create a stub structure.
static Vector< T > solve(const SparseMatrix< T, impl > &A, const Vector< T > &b)
Solve the linear system.
[Definition of the system]
static float val()
Eta is constant one.
void const_field
indicate it is a constant field
[Definition of the system]
It ancapsulate the minus operation.
Definition sum.hpp:142
It model an expression expr1 * expr2.
Definition mul.hpp:120
It model an expression expr1 + ... exprn.
Definition sum.hpp:93