OpenFPM_pdata  4.1.0 Project that contain the implementation of distributed structures
main_eigen.cpp
1
40
42 #include "Grid/grid_dist_id.hpp"
43 #include "Matrix/SparseMatrix.hpp"
44 #include "Vector/Vector.hpp"
45 #include "FiniteDifference/FDScheme.hpp"
46 #include "FiniteDifference/util/common.hpp"
47 #include "FiniteDifference/eq.hpp"
48 #include "Solvers/umfpack_solver.hpp"
49 #include "Solvers/petsc_solver.hpp"
50
51 struct lid_nn
52 {
53  // dimensionaly of the equation (2D problem 3D problem ...)
54  static const unsigned int dims = 2;
55
56  // number of fields in the system v_x, v_y, P so a total of 3
57  static const unsigned int nvar = 3;
58
59  // boundary conditions PERIODIC OR NON_PERIODIC
60  static const bool boundary[];
61
62  // type of space float, double, ...
63  typedef float stype;
64
65  // type of base grid, it is the distributed grid that will store the result
66  // Note the first property is a 2D vector (velocity), the second is a scalar (Pressure)
68
69  // type of SparseMatrix, for the linear system, this parameter is bounded by the solver
70  // that you are using, in case of umfpack using <double,int> it is the only possible choice
71  // By default SparseMatrix is EIGEN based
72  typedef SparseMatrix<double,int> SparseMatrix_type;
73
74  // type of Vector for the linear system, this parameter is bounded by the solver
75  // that you are using, in case of umfpack using <double> it is the only possible choice
76  typedef Vector<double> Vector_type;
77
78  // Define that the underline grid where we discretize the system of equation is staggered
79  static const int grid_type = STAGGERED_GRID;
80 };
81
82 const bool lid_nn::boundary[] = {NON_PERIODIC,NON_PERIODIC};
83
85
86
117
119 // Constant Field
120 struct eta
121 {
122  typedef void const_field;
123
124  static float val() {return 1.0;}
125 };
126
127 // Convenient constants
128 constexpr unsigned int v[] = {0,1};
129 constexpr unsigned int P = 2;
130 constexpr unsigned int ic = 2;
131 constexpr int x = 0;
132 constexpr int y = 1;
133
134 // Create field that we have v_x, v_y, P
135 typedef Field<v[x],lid_nn> v_x; // Definition 1 v_x
136 typedef Field<v[y],lid_nn> v_y; // Definition 2 v_y
137 typedef Field<P,lid_nn> Prs; // Definition 3 Prs
138
139 // Eq1 V_x
140
141 typedef mul<eta,Lap<v_x,lid_nn>,lid_nn> eta_lap_vx; // Step 1
142 typedef D<x,Prs,lid_nn> p_x; // Step 2
143 typedef minus<p_x,lid_nn> m_p_x; // Step 3
144 typedef sum<eta_lap_vx,m_p_x,lid_nn> vx_eq; // Step 4
145
146 // Eq2 V_y
147
149 typedef D<y,Prs,lid_nn> p_y;
150 typedef minus<p_y,lid_nn> m_p_y;
152
153 // Eq3 Incompressibility
154
155 typedef D<x,v_x,lid_nn,FORWARD> dx_vx; // Step 5
156 typedef D<y,v_y,lid_nn,FORWARD> dy_vy; // Step 6
157 typedef sum<dx_vx,dy_vy,lid_nn> ic_eq; // Step 7
158
160
193
195 // Equation for boundary conditions
196
197 // Directional Avg
198 typedef Avg<x,v_y,lid_nn> avg_vy;
199 typedef Avg<y,v_x,lid_nn> avg_vx;
200
203
204 // Usefull constants (as MACRO)
205 #define EQ_1 0
206 #define EQ_2 1
207 #define EQ_3 2
208
210
211 #include "Vector/vector_dist.hpp"
212 #include "data_type/aggregate.hpp"
213
214 int main(int argc, char* argv[])
215 {
235
237  // Initialize
238  openfpm_init(&argc,&argv);
239
240  // velocity in the grid is the property 0, pressure is the property 1
241  constexpr int velocity = 0;
242  constexpr int pressure = 1;
243
244  // Domain, a rectangle
245  Box<2,float> domain({0.0,0.0},{3.0,1.0});
246
247  // Ghost (Not important in this case but required)
248  Ghost<2,float> g(0.01);
249
250  // Grid points on x=256 and y=64
251  long int sz[] = {96,32};
252  size_t szu[2];
253  szu[0] = (size_t)sz[0];
254  szu[1] = (size_t)sz[1];
255
256  // We need one more point on the left and down part of the domain
257  // This is given by the boundary conditions that we impose.
258  //
260
262
274
277
279
297
299  // It is the maximum extension of the stencil (order 2 laplacian stencil has extension 1)
300  Ghost<2,long int> stencil_max(1);
301
302  // Finite difference scheme
303  FDScheme<lid_nn> fd(pd, stencil_max, domain, g_dist);
304
306
330
332  fd.impose(ic_eq(),0.0, EQ_3, {0,0},{sz[0]-2,sz[1]-2},true);
333  fd.impose(Prs(), 0.0, EQ_3, {0,0},{0,0});
334
335  // Here we impose the Eq1 and Eq2
336  fd.impose(vx_eq(),0.0, EQ_1, {1,0},{sz[0]-2,sz[1]-2});
337  fd.impose(vy_eq(),0.0, EQ_2, {0,1},{sz[0]-2,sz[1]-2});
338
339  // v_x and v_y
340  // Imposing B1
341  fd.impose(v_x(),0.0, EQ_1, {0,0},{0,sz[1]-2});
342  fd.impose(avg_vy_f(),0.0, EQ_2 , {-1,0},{-1,sz[1]-1});
343  // Imposing B2
344  fd.impose(v_x(),0.0, EQ_1, {sz[0]-1,0},{sz[0]-1,sz[1]-2});
345  fd.impose(avg_vy(),1.0, EQ_2, {sz[0]-1,0},{sz[0]-1,sz[1]-1});
346
347  // Imposing B3
348  fd.impose(avg_vx_f(),0.0, EQ_1, {0,-1},{sz[0]-1,-1});
349  fd.impose(v_y(), 0.0, EQ_2, {0,0},{sz[0]-2,0});
350  // Imposing B4
351  fd.impose(avg_vx(),0.0, EQ_1, {0,sz[1]-1},{sz[0]-1,sz[1]-1});
352  fd.impose(v_y(), 0.0, EQ_2, {0,sz[1]-1},{sz[0]-2,sz[1]-1});
353
354  // When we pad the grid, there are points of the grid that are not
355  // touched by the previous condition. Mathematically this lead
356  // to have too many variables for the conditions that we are imposing.
357  // Here we are imposing variables that we do not touch to zero
358  //
359
361  fd.impose(Prs(), 0.0, EQ_3, {-1,-1},{sz[0]-1,-1});
362  fd.impose(Prs(), 0.0, EQ_3, {-1,sz[1]-1},{sz[0]-1,sz[1]-1});
363  fd.impose(Prs(), 0.0, EQ_3, {-1,0},{-1,sz[1]-2});
364  fd.impose(Prs(), 0.0, EQ_3, {sz[0]-1,0},{sz[0]-1,sz[1]-2});
365
367  fd.impose(v_x(), 0.0, EQ_1, {-1,-1},{-1,sz[1]-1});
368  fd.impose(v_y(), 0.0, EQ_2, {-1,-1},{sz[0]-1,-1});
369
371
385
387  // Create an UMFPACK solver
388  umfpack_solver<double> solver;
389
390  // Give to the solver A and b, return x, the solution
391  auto x = solver.solve(fd.getA(),fd.getB());
392
394
406
408  fd.template copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1},g_dist);
409
410  g_dist.write("lid_driven_cavity_p_umfpack");
411
413
414
426
428  openfpm_finalize();
429
431
432
441 }
Derivative second order on h (spacing)
Definition: Derivative.hpp:28
[Definition of the system]
It ancapsulate the minus operation.
Definition: sum.hpp:141
Sparse Matrix implementation stub object when OpenFPM is compiled with no linear algebra support.
Definition: Vector.hpp:39
Definition: eq.hpp:82
Class that contain Padding information on each direction positive and Negative direction.
Definition: Ghost.hpp:126
stub when library compiled without eigen
Definition: Ghost.hpp:39
It model an expression expr1 * expr2.
Definition: mul.hpp:119
This class decompose a space into sub-sub-domains and distribute them across processors.
Sparse Matrix implementation.
This is a distributed grid.
Average.
Definition: Average.hpp:26
[Definition of the system]
This class represent an N-dimensional box.
Definition: Box.hpp:60
static Vector< double > solve(SparseMatrix< double, id_type, impl > &A, const Vector< double > &b, size_t opt=UMFPACK_NONE)
stub solve
It model an expression expr1 + ... exprn.
Definition: sum.hpp:92
Finite Differences.
Definition: FDScheme.hpp:126
Test structure used for several test.
Definition: Point_test.hpp:105