OpenFPM_pdata  4.1.0 Project that contain the implementation of distributed structures main_petsc.cpp
1
31
33 #include "config.h"
34
35 #ifdef HAVE_PETSC
36
37 #include "Grid/grid_dist_id.hpp"
38 #include "Matrix/SparseMatrix.hpp"
39 #include "Vector/Vector.hpp"
40 #include "FiniteDifference/FDScheme.hpp"
41 #include "FiniteDifference/util/common.hpp"
42 #include "FiniteDifference/eq.hpp"
43 #include "Solvers/petsc_solver.hpp"
44 #include "Solvers/petsc_solver.hpp"
45
46 struct lid_nn
47 {
48  // dimensionaly of the equation (2D problem 3D problem ...)
49  static const unsigned int dims = 2;
50
51  // number of fields in the system v_x, v_y, P so a total of 3
52  static const unsigned int nvar = 3;
53
54  // boundary conditions PERIODIC OR NON_PERIODIC
55  static const bool boundary[];
56
57  // type of space float, double, ...
58  typedef float stype;
59
60  // type of base grid, it is the distributed grid that will store the result
61  // Note the first property is a 2D vector (velocity), the second is a scalar (Pressure)
63
64  // type of SparseMatrix, for the linear system, this parameter is bounded by the solver
65  // that you are using, in case of umfpack using <double,int> it is the only possible choice
66  // Here we choose PETSC implementation
67  typedef SparseMatrix<double,int,PETSC_BASE> SparseMatrix_type;
68
69  // type of Vector for the linear system, this parameter is bounded by the solver
70  // that you are using
71  typedef Vector<double,PETSC_BASE> Vector_type;
72
73  // Define that the underline grid where we discretize the system of equation is staggered
74  static const int grid_type = STAGGERED_GRID;
75 };
76
77 const bool lid_nn::boundary[] = {NON_PERIODIC,NON_PERIODIC};
78
80
81
112
114 // Constant Field
115 struct eta
116 {
117  typedef void const_field;
118
119  static float val() {return 1.0;}
120 };
121
122 // Convenient constants
123 constexpr unsigned int v[] = {0,1};
124 constexpr unsigned int P = 2;
125 constexpr unsigned int ic = 2;
126 constexpr int x = 0;
127 constexpr int y = 1;
128
129 // Create field that we have v_x, v_y, P
130 typedef Field<v[x],lid_nn> v_x; // Definition 1 v_x
131 typedef Field<v[y],lid_nn> v_y; // Definition 2 v_y
132 typedef Field<P,lid_nn> Prs; // Definition 3 Prs
133
134 // Eq1 V_x
135
136 typedef mul<eta,Lap<v_x,lid_nn>,lid_nn> eta_lap_vx; // Step 1
137 typedef D<x,Prs,lid_nn> p_x; // Step 2
138 typedef minus<p_x,lid_nn> m_p_x; // Step 3
139 typedef sum<eta_lap_vx,m_p_x,lid_nn> vx_eq; // Step 4
140
141 // Eq2 V_y
142
144 typedef D<y,Prs,lid_nn> p_y;
145 typedef minus<p_y,lid_nn> m_p_y;
147
148 // Eq3 Incompressibility
149
150 typedef D<x,v_x,lid_nn,FORWARD> dx_vx; // Step 5
151 typedef D<y,v_y,lid_nn,FORWARD> dy_vy; // Step 6
152 typedef sum<dx_vx,dy_vy,lid_nn> ic_eq; // Step 7
153
155
188
190 // Equation for boundary conditions
191
192 // Directional Avg
193 typedef Avg<x,v_y,lid_nn> avg_vy;
194 typedef Avg<y,v_x,lid_nn> avg_vx;
195
198
199 // Usefull constants (as MACRO)
200 #define EQ_1 0
201 #define EQ_2 1
202 #define EQ_3 2
203
205
206 #include "Vector/vector_dist.hpp"
207 #include "data_type/aggregate.hpp"
208
209 int main(int argc, char* argv[])
210 {
230
232  // Initialize
233  openfpm_init(&argc,&argv);
234
235  // velocity in the grid is the property 0, pressure is the property 1
236  constexpr int velocity = 0;
237  constexpr int pressure = 1;
238
239  // Domain, a rectangle
240  Box<2,float> domain({0.0,0.0},{12.0,4.0});
241
242  // Ghost (Not important in this case but required)
243  Ghost<2,float> g(0.01);
244
245  // Grid points on x=96 and y=32
246  long int sz[] = {96,32};
247  size_t szu;
248  szu = (size_t)sz;
249  szu = (size_t)sz;
250
251  // We need one more point on the left and down part of the domain
252  // This is given by the boundary conditions that we impose.
253  //
255
257
269
272
274
292
294  // It is the maximum extension of the stencil (order 2 laplacian stencil has extension 1)
295  Ghost<2,long int> stencil_max(1);
296
297  // Finite difference scheme
298  FDScheme<lid_nn> fd(pd, stencil_max, domain, g_dist);
299
301
325
327  fd.impose(ic_eq(),0.0, EQ_3, {0,0},{sz-2,sz-2},true);
328  fd.impose(Prs(), 0.0, EQ_3, {0,0},{0,0});
329
330  // Here we impose the Eq1 and Eq2
331  fd.impose(vx_eq(),0.0, EQ_1, {1,0},{sz-2,sz-2});
332  fd.impose(vy_eq(),0.0, EQ_2, {0,1},{sz-2,sz-2});
333
334  // v_x and v_y
335  // Imposing B1
336  fd.impose(v_x(),0.0, EQ_1, {0,0},{0,sz-2});
337  fd.impose(avg_vy_f(),0.0, EQ_2 , {-1,0},{-1,sz-1});
338  // Imposing B2
339  fd.impose(v_x(),0.0, EQ_1, {sz-1,0},{sz-1,sz-2});
340  fd.impose(avg_vy(),1.0, EQ_2, {sz-1,0},{sz-1,sz-1});
341
342  // Imposing B3
343  fd.impose(avg_vx_f(),0.0, EQ_1, {0,-1},{sz-1,-1});
344  fd.impose(v_y(), 0.0, EQ_2, {0,0},{sz-2,0});
345  // Imposing B4
346  fd.impose(avg_vx(),0.0, EQ_1, {0,sz-1},{sz-1,sz-1});
347  fd.impose(v_y(), 0.0, EQ_2, {0,sz-1},{sz-2,sz-1});
348
349  // When we pad the grid, there are points of the grid that are not
350  // touched by the previous condition. Mathematically this lead
351  // to have too many variables for the conditions that we are imposing.
352  // Here we are imposing variables that we do not touch to zero
353  //
354
356  fd.impose(Prs(), 0.0, EQ_3, {-1,-1},{sz-1,-1});
357  fd.impose(Prs(), 0.0, EQ_3, {-1,sz-1},{sz-1,sz-1});
358  fd.impose(Prs(), 0.0, EQ_3, {-1,0},{-1,sz-2});
359  fd.impose(Prs(), 0.0, EQ_3, {sz-1,0},{sz-1,sz-2});
360
362  fd.impose(v_x(), 0.0, EQ_1, {-1,-1},{-1,sz-1});
363  fd.impose(v_y(), 0.0, EQ_2, {-1,-1},{sz-1,-1});
364
366
385
387  // Create a PETSC solver
388  petsc_solver<double> solver;
389
390  // Set the maxumum nunber if iterations
391  solver.setMaxIter(1000);
392
393  solver.setRestart(200);
394
395  // Give to the solver A and b, return x, the solution
396  auto x = solver.try_solve(fd.getA(),fd.getB());
397
399
411
413  fd.template copy<velocity,pressure>(x,{0,0},{sz-1,sz-1},g_dist);
414
415  g_dist.write("lid_driven_cavity_p_petsc");
416
418
419
431
433  openfpm_finalize();
434
436
437
446 }
447
448
449 #else
450
451 int main(int argc, char* argv[])
452 {
453  return 0;
454 }
455
456 #endif
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
void setMaxIter(PetscInt n)
Set the maximum number of iteration for Krylov solvers.
Class that contain Padding information on each direction positive and Negative direction.
Definition: Ghost.hpp:126
Definition: Ghost.hpp:39
Vector< double, PETSC_BASE > try_solve(SparseMatrix< double, int, PETSC_BASE > &A, const Vector< double, PETSC_BASE > &b)
Try to solve the system using all the solvers and generate a report.
void setRestart(PetscInt n)
For GMRES based method, the number of Krylov directions to orthogonalize against.
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.
This class is able to do Matrix inversion in parallel with PETSC solvers.
Average.
Definition: Average.hpp:26
[Definition of the system]
This class represent an N-dimensional box.
Definition: Box.hpp:60
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